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

[Matthew Knepley]

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/20150909/047466bc/attachment-0001.html>