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

[Justin Chang]

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?

Thanks,
Justin
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150909/f8094b98/attachment.html>