<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Thu, Sep 10, 2015 at 5:49 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">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.<div><br></div><div>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:</div></div></blockquote><div><br></div><div>I do not see the advantage of a variational form here. There is no "energy" for the problem. Colocation seems easier.<br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div>#==========================</div><div>mesh = ...</div><div>Q = FunctionSpace(mesh, "CG", 1)</div><div>W = Q*Q*Q</div><div>dc_A,dc_B,dc_C = TrialFunctions(W)</div><div>q_A,q_B,q_C = TestFunctions(W)</div><div>c = Function(W)</div><div>c_A,c_B,c_C = c.split()</div><div><br></div><div>...</div><div><br></div><div>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</div><div>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</div><div><br></div><div>...</div><div><br></div><div># Solve for u0_A a u0_B</div><div><br></div><div>....</div><div><br></div><div>solve(F == 0, c, J=J)</div><div>#===========================</div></div><div><br></div><div>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?</div><div><br></div><div>Or is this even the right approach?</div></div></blockquote><div><br></div><div>I would just do simple FD (colocation) for this part of the problem. Then you can see everything easily. Use</div><div>the FD Jacobian at first. You should converge if your guess is good enough. Then put in an analytic Jacobian.</div><div><br></div><div>  Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Thanks,</div><div>Justin</div><div><br></div></div><div class="HOEnZb"><div class="h5"><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Sep 9, 2015 at 4:41 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><div>On Wed, Sep 9, 2015 at 4:34 PM, Justin Chang <span dir="ltr"><<a href="mailto:jychang48@gmail.com" target="_blank">jychang48@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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></div></blockquote><div><br></div></div></div><div>1) I had the same thing in DFT where I had to solve a local equation for the screening length at each node</div><div><br></div><div>2) I would use a SNES</div><div><br></div><div>3) If you think there is variability in the solve, then do each individually. Otherwise,</div><div>    you can easily do a Newton with a diagonal Jacobian.</div><div><br></div><div>4) I would eliminate c_C since its so easy</div><div><br></div><div>5) These look quadratic. Can't you solve this analytically?</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><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><span><font color="#888888"><br></font></span></div><span><font color="#888888">
</font></span></blockquote></div><span><font color="#888888"><br><br clear="all"><div><br></div>-- <br><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</font></span></div></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>