Never mind I think I figured it out. Will let you know if I think something is wrong. Thanks anyway
On Wed, May 6, 2015 at 8:06 AM, Matthew Knepley <[email protected]> wrote: > On Wed, May 6, 2015 at 2:07 AM, Justin Chang <[email protected]> wrote: > >> Matt, >> >> 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. >> >> 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. >> >> 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. >> > > For some reason, I thought that was taken care of. Please resend. Sorry > about that. > > Thanks, > > Matt > > >> Thanks, >> Justin >> >> On Tue, May 5, 2015 at 3:34 PM, Matthew Knepley <[email protected]> >> wrote: >> >>> On Sat, Apr 4, 2015 at 11:54 AM, Justin Chang <[email protected]> wrote: >>> >>>> Hi everyone, >>>> >>>> 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: >>>> >>>> du/dt + v * grad[c] - div[grad[c]] = f >>>> >>>> Where v is the velocity (obtained from the auxiliary term a_x[]). For >>>> the residual, would it be of the following form: >>>> >>> >>> 1) I would think you would use a[] instead. What is your velocity the >>> gradient of? >>> >>> >>>> 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[]) { >>>> >>>> PetscInt d; >>>> f[0] = u_t[0]; >>>> for (d = 0; d < spatialDim; ++d) f0[0] += a_x[d]*u_x[d]; >>>> >>>> } >>>> >>>> What about the jacobian? My guess would be to use g1(...) but what >>>> would the inside of this function be? >>>> >>> >>> 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} >>> for 2D, so g1 would have two terms, >>> >>> g1[0] = a[0]; >>> g1[1] = a[1]; >>> >>> Thanks, >>> >>> Matt >>> >>> >>>> Thanks, >>>> >>>> >>>> -- >>>> Justin Chang >>>> PhD Candidate, Civil Engineering - Computational Sciences >>>> University of Houston, Department of Civil and Environmental Engineering >>>> Houston, TX 77004 >>>> (512) 963-3262 >>>> >>> >>> >>> >>> -- >>> 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 >
