<div dir="ltr"><div><div><div><div>Matt,<br><br></div>Yes I would actually use a[] in this case. Back then, I think I was trying to take the derivative of the auxiliary which i had declared to be piece-wise constant.<br><br>The auxiliary a[] is the gradient of pressure (times permeability over viscosity) obtained from solving the Darcy equation as a Laplacian. This part is solved separately, projected as cell-wise velocity, and inputted as an auxiliary for the advection-diffusion problem.<br><br></div>Speaking of which, have you had a chance to look into DMPlexProjectFieldLocal(...) not projecting dirichlet BC constraints? I realize that the "small example" I gave you wasn't small at all, so I can send you a smaller and simpler one if needed.<br><br></div>Thanks,<br></div>Justin<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, May 5, 2015 at 3:34 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"><span class="">On Sat, Apr 4, 2015 at 11:54 AM, Justin Chang <span dir="ltr"><<a href="mailto:jchang27@uh.edu" target="_blank">jchang27@uh.edu</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><div><div><div><div>Hi everyone,<br><br></div>I want to include advection into my diffusion FEM code (I am assuming a small Peclet number so stability isn't an issue ATM). That is I want to incorporate the second term as a pointwise function:<br><br></div>du/dt + v * grad[c] - div[grad[c]] = f<br><br></div>Where v is the velocity (obtained from the auxiliary term a_x[]). For the residual, would it be of the following form:<br></div></div></div></blockquote><div><br></div></span><div>1) I would think you would use a[] instead. What is your velocity the gradient of?</div><span class=""><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>f0_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]) {<br><br></div><div>  PetscInt d;<br></div><div>  f[0] = u_t[0];<br></div>  for (d = 0; d < spatialDim; ++d) f0[0] += a_x[d]*u_x[d];<br><div><br>}<br><br></div><div>What about the jacobian? My guess would be to use g1(...) but what would the inside of this function be?<br></div></div></blockquote><div><br></div></span><div>Yes it would be g1. The indices for the output are f,g,dg. I am guessing that c is a scalar, so f = {0}, g = {0}, dg = {0, 1}</div><div>for 2D, so g1 would have two terms,</div><div><br></div><div>  g1[0] = a[0];</div><div>  g1[1] = a[1];</div><div><br></div><div>  Thanks,</div><div><br></div><div>    Matt</div><span class=""><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>Thanks,<span><font color="#888888"><br></font></span></div><span><font color="#888888"><br clear="all"><br>-- <br><div><div dir="ltr"><div><div><div>Justin Chang<br></div>PhD Candidate, Civil Engineering - Computational Sciences<br></div>University of Houston, Department of Civil and Environmental Engineering<br></div>Houston, TX 77004<br><a href="tel:%28512%29%20963-3262" value="+15129633262" target="_blank">(512) 963-3262</a><br></div></div></font></span></div>
</blockquote></span></div><span class="HOEnZb"><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>