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

[Matthew Knepley]

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/47e6039e/attachment-0001.html>