Hi all,

Hmm strangely that did the trick. Thanks!

As for the "accuracy" of these concentrations going into the reactions, that's something I will need to test thoroughly since I am doing operator splitting of A and D and not solving an actual AD equation. But it is known that anisotropy and advection violates the discrete maximum principles (i.e., negative concentrations) so one of our ways to address this issue is through optimization.

Justin

On Fri, Sep 11, 2015 at 3:55 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,

> On 11 Sep 2015, at 00:56, Justin Chang <jychang48@gmail.com> wrote:
>
> Hi everyone,
>
> Attached is a tentative working code I have for an advective-diffusive-reactive system using operator splitting. I tried to express the system of equations as a mixed finite element form (not sure if this is the right approach).
>
> When I run the attached code, i get an error saying "Nonlinear solver failed to converge after 0 nonlinear iterations. Reason: DIVERGED_LINE_SEARCH".
> Can you guys have a look and see if this is even a logical approach? The parts in question begin at line 133, and the solver is at line 217 (if you comment out line 217, the advection-diffusion part runs well).
>
> If this doesn't work, I will explore the pyop2 option.

You've been bitten by the most subtle part of the fenics interface!  You're not the first.  https://github.com/firedrakeproject/firedrake/issues/407 and https://bitbucket.org/fenics-project/dolfin/issues/194/split-u-versus-usplit

Specifically, there's an unfortunate difference in the behaviour of the results of:

foo = Function(MixedFunctionSpace)

a, b = foo.split()

as opposed to:

a, b = split(foo)

The former is useable in outputs (and for pointwise assignments), however, they don't function as fully-fledged "coefficients" in forms.

The latter is /not/ useable in outputs (and pointless assignments).  But does function correctly in forms.

I was initially confused as to why you'd hand-coded the Jacobian for the species reaction at all, so I tried taking it out (in case it was wrong), and got this error:

/Users/lmitche1/Documents/work/fenics/ufl/ufl/log.pyc in error(self, *message)
    149         "Write error message and raise an exception."
    150         self._log.error(*message)
--> 151         raise self._exception_type(self._format_raw(*message))
    152
    153     def begin(self, *message):

UFLException: Expecting all integrals in a form to share geometric dimension, got ().

Which you probably saw as well.

This occurs because you've used the first split form above, and so UFL doesn't know about the split coefficients in the form and things it can't differentiate the residual with respect to the full state c.


If I go ahead and change the version of split I use, everything works nicely (last least I don't get convergence errors).

W_R = Q*Q*Q
dc_A,dc_B,dc_C = TrialFunctions(W_R)
q_A,q_B,q_C = TestFunctions(W_R)
c = Function(W_R)
# Use UFL's "split" here so that the residual knows about where
# the components come from
c_A,c_B,c_C = split(c)

# Residual and Jacobian form
F = (q_A*(c_A + c_C - u0_A) + q_B*(c_B + c_C - u0_B) + q_C*(c_A*c_B - c_C))*dx
problem_R = NonlinearVariationalProblem(F,c)
solver_R = NonlinearVariationalSolver(problem_R,solver_parameters=R_parameters,options_prefix="R_")
# Now go back to the firedrake "split" so that we can do pointless
# assignments later on.
c_A,c_B,c_C = c.split()


Cheers,

Lawrence

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake