[petsc-users] Expressing this nonlinear system in PETSc and/or petsc4py

[Justin Chang]

Matt,

So it turns out that the issue I had was using c.split() as opposed to
split(c). The firedrake folks helped me resolve this (for now).

Thanks,
Justin

On Fri, Sep 11, 2015 at 5:45 AM, Matthew Knepley <knepley at gmail.com> wrote:

> On Thu, Sep 10, 2015 at 5:49 PM, Justin Chang <jychang48 at gmail.com> wrote:
>
>> The example I have happens to be a simple system (e.g., Na^(+) + Cl^(-)
>> <=> NaCl), I could potentially solve a 7 component system with 4 primary
>> and 3 secondary species which in some cases cannot be solved analytically
>> as easily.
>>
>> Anyway, I was wondering if I could express the example system I had as a
>> mixed weak/variational form. That is, multiply each equation by a test
>> function and solve this system as if it were a finite element problem. I
>> had the following weak/variational form expressed in the FEniCS/Firedrake
>> form:
>>
>
> I do not see the advantage of a variational form here. There is no
> "energy" for the problem. Colocation seems easier.
>
>
>> #==========================
>> mesh = ...
>> Q = FunctionSpace(mesh, "CG", 1)
>> W = Q*Q*Q
>> dc_A,dc_B,dc_C = TrialFunctions(W)
>> q_A,q_B,q_C = TestFunctions(W)
>> c = Function(W)
>> c_A,c_B,c_C = c.split()
>>
>> ...
>>
>> F = (q_A*(u0_A - c_A - c_C) + q_B*(u0_B - c_B - c_C) + q_C*(c_C -
>> c_A*c_B))*dx
>> J = (-q_A*(dc_A + dc_C) + q_B*(dc_B + dc_C) + q_C*(dc_C - c_B*dc_A -
>> c_A*dc_B))*dx
>>
>> ...
>>
>> # Solve for u0_A a u0_B
>>
>> ....
>>
>> solve(F == 0, c, J=J)
>> #===========================
>>
>> Where u0_A and u0_B correspond to my psi_A/B from a different function
>> space. However, when I solve the system, I get a DIVERGED_LINE_SEARCH. From
>> what I have read, it seems this error mostly arises because of an incorrect
>> Jacobian. But I think my Jacobian J looks correct, is it not?
>>
>> Or is this even the right approach?
>>
>
> I would just do simple FD (colocation) for this part of the problem. Then
> you can see everything easily. Use
> the FD Jacobian at first. You should converge if your guess is good
> enough. Then put in an analytic Jacobian.
>
>   Matt
>
>
>> Thanks,
>> Justin
>>
>>
>> On Wed, Sep 9, 2015 at 4:41 PM, Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Wed, Sep 9, 2015 at 4:34 PM, Justin Chang <jychang48 at gmail.com>
>>> wrote:
>>>
>>>> Hi everyone,
>>>>
>>>> I need to solve this system of nonlinear geochemical reactions:
>>>>
>>>> psi_A = c_A + c_C
>>>> psi_B = c_B + c_C
>>>> c_C = k*c_A*c_B
>>>>
>>>> where psi_A and psi_B are sub components (i.e., a mixed system) of my
>>>> mesh (this was done using firedrake) and k is a constant scalar. These
>>>> quantities are known a prior, and psi_A/B was obtained from the
>>>> advection-diffusion equation (using SUPG). Given the two-field formulation
>>>> of psi_A/B, I need to obtain a three-field formulation of c_A/B/C using the
>>>> above system of equations. It should be noted that the above system of
>>>> equations are node-independent. That is, c_A/B/C of one node does not care
>>>> what psi_A/B of other nodes may be.
>>>>
>>>> What’s the best strategy to go about solving this? With SciPy, i did
>>>> something like this:
>>>>
>>>> #====================
>>>> from scipy.optimize import fsolve
>>>> import math
>>>>
>>>> # advection-diffusion for psi_A and psi_B
>>>>
>>>> # Nonlinear function for geochemical reactions
>>>> def equations(p,psi_A,psi_B,k_1):
>>>>  c_A,c_B,c_C = p
>>>>  return (c_A+c_C-psi_A,c_B+c_C-psi_B,c_C-k_1*c_A*c_B)
>>>>
>>>> # Initialize
>>>> c_A_vec = np.zeros(len(psi_A.vector()))
>>>> c_B_vec = np.zeros(len(psi_A.vector()))
>>>> c_C_vec = np.zeros(len(psi_A.vector()))
>>>> k_1 = 1.0
>>>>
>>>> for i in range(len(psi_A.vector())):
>>>>  c_A,c_B,c_C =
>>>> fsolve(equations,(0.0,0.0,0.0),args=(psi_A.vector()[i],psi_B.vector()[i],k_1))
>>>>  c_A_vec[i] = c_A
>>>>  c_B_vec[i] = c_B
>>>>  c_C=vec[i] = c_C
>>>>
>>>> c_A = Function(Q)
>>>> c_B = Function(Q)
>>>> c_C = Function(Q)
>>>> c_A.vector()[:]=c_A_vec
>>>> c_B.vector()[:]=c_B_vec
>>>> c_C.vector()[:]=c_C_vec
>>>> #======================
>>>>
>>>> The above is my temporary work-around to this issue. Basically, what I
>>>> did was I solved the equations at each node. But is there a PETSc or
>>>> petsc4py way to do this, specifically solving the above equations globally
>>>> instead of at each individual node?
>>>>
>>>
>>> 1) I had the same thing in DFT where I had to solve a local equation for
>>> the screening length at each node
>>>
>>> 2) I would use a SNES
>>>
>>> 3) If you think there is variability in the solve, then do each
>>> individually. Otherwise,
>>>     you can easily do a Newton with a diagonal Jacobian.
>>>
>>> 4) I would eliminate c_C since its so easy
>>>
>>> 5) These look quadratic. Can't you solve this analytically?
>>>
>>>   Thanks,
>>>
>>>      Matt
>>>
>>>
>>>>
>>>> Thanks,
>>>> Justin
>>>>
>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150911/d1f27aca/attachment.html>