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 <[email protected]> wrote: > On Wed, Sep 9, 2015 at 4:34 PM, Justin Chang <[email protected]> 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 >
