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

[Justin Chang]

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:

#==========================
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?

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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150910/c846ab88/attachment.html>