Matt, So it turns out that the issue I had was using c.split() as opposed to split(c). The firedrake folks helped me resolve this (for now).
Thanks, Justin On Fri, Sep 11, 2015 at 5:45 AM, Matthew Knepley <[email protected]> wrote: > On Thu, Sep 10, 2015 at 5:49 PM, Justin Chang <[email protected]> wrote: > >> 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: >> > > I do not see the advantage of a variational form here. There is no > "energy" for the problem. Colocation seems easier. > > >> #========================== >> 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? >> > > I would just do simple FD (colocation) for this part of the problem. Then > you can see everything easily. Use > the FD Jacobian at first. You should converge if your guess is good > enough. Then put in an analytic Jacobian. > > Matt > > >> 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 >>> >> >> > > > -- > 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 >
