<div dir="ltr"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">Hi everyone,</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">I need to solve this system of nonlinear geochemical reactions:</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">psi_A = c_A + c_C</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">psi_B = c_B + c_C</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_C = k*c_A*c_B</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">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.</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">What’s the best strategy to go about solving this? With SciPy, i did something like this:</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">#====================</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">from scipy.optimize import fsolve</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">import math</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"># advection-diffusion for psi_A and psi_B</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"># Nonlinear function for geochemical reactions</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">def equations(p,psi_A,psi_B,k_1):</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"> c_A,c_B,c_C = p </span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"> return (c_A+c_C-psi_A,c_B+c_C-psi_B,c_C-k_1*c_A*c_B)</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"># Initialize</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_A_vec = np.zeros(len(psi_A.vector()))</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_B_vec = np.zeros(len(psi_A.vector()))</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_C_vec = np.zeros(len(psi_A.vector()))</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">k_1 = 1.0 </span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">for i in range(len(psi_A.vector())):</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"> 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))</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"> c_A_vec[i] = c_A</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"> c_B_vec[i] = c_B</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"> c_C=vec[i] = c_C </span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_A = Function(Q)</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_B = Function(Q)</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_C = Function(Q)</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_A.vector()[:]=c_A_vec</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_B.vector()[:]=c_B_vec</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">c_C.vector()[:]=c_C_vec</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">#======================</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">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?</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">Thanks,</span><br style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px"><span style="color:rgb(0,0,0);font-family:Helvetica;font-size:12px">Justin</span><br></div>