[petsc-users] DMPlexProjectField not including BC constraints?

[Matthew Knepley]

On Wed, Apr 22, 2015 at 9:57 PM, Justin Chang <jchang27 at uh.edu> wrote:

> Hello,
>
> I am using DMPlexProjectField(...) to postprocess my final solution u into
> another vector sol. I have the following pointwise functions:
>
> /* x-component of gradient */
> void velx(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 v[])
> {
>   v[0] = -a[0]*mu*u_x[0];
> }
>
> /* y-component of gradient */
> void vely(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 v[])
> {
>   v[0] = -a[0]*mu*u_x[1];
> }
>
> /* z-component of gradient */
> void velz(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 v[])
> {
>   v[0] = -a[0]*mu*u_x[2];
> }
>
> Where a[0] is the cell-wise permeability and mu is the (inverse of)
> viscosity. I am solving the diffusion equation (P1 elements). At the end of
> my simulation I have these lines:
>
> ierr = VecDuplicate(user.u,&user.sol);CHKERRQ(ierr);
> ierr = PetscObjectSetName((PetscObject) user.sol,
> "velocity");CHKERRQ(ierr);
> void (*vX[1])(const PetscScalar *, const PetscScalar *, const PetscScalar
> *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const
> PetscReal *, PetscScalar *) = {velx};
> ierr =
> DMPlexProjectField(dm,user.u,vX,INSERT_ALL_VALUES,user.sol);CHKERRQ(ierr);
> ierr = VecViewFromOptions(user.sol,NULL,"-vec_view_x");CHKERRQ(ierr);
>
> void (*vY[1])(const PetscScalar *, const PetscScalar *, const PetscScalar
> *, const PetscScalar *, const PetscScalar *, const PetscScalar *, const
> PetscReal *, PetscScalar *) = {vely};
> ierr =
> DMPlexProjectField(dm,user.u,vY,INSERT_ALL_VALUES,user.sol);CHKERRQ(ierr);
> ierr = VecViewFromOptions(user.sol,NULL,"-vec_view_y");CHKERRQ(ierr);
>
> if (spatialDim == 3) {
>   void (*vZ[1])(const PetscScalar *, const PetscScalar *, const
> PetscScalar *, const PetscScalar *, const PetscScalar *, const PetscScalar
> *, const PetscReal *, PetscScalar *) = {velz};
>   ierr =
> DMPlexProjectField(dm,user.u,vZ,INSERT_ALL_VALUES,user.sol);CHKERRQ(ierr);
>   ierr = VecViewFromOptions(user.sol,NULL,"-vec_view_z");CHKERRQ(ierr);
> }
> ierr = VecDestroy(&user.sol);CHKERRQ(ierr);
>
> But when i look at the results, the dirichlet constrained values of u were
> projected into the respective locations in sol. I thought putting
> INSERT_ALL_VALUES would insert the pointwise functions into all locations
> including the constrained ones. How can I project the above pointwise
> functions into the constrained fields? Or what am I potentially missing
> here?
>

There are no Dirichlet constrained values in the global vector user.sol.
When you call VecView() it gets projected to
an enlarged vector including the BC and output. What you want is to project
that vector into a space with no constraints,
so you really need another DM/Section combination.

I think this will work:

  1) Get a local vector from dm

  2) Use DMPlexProjectFieldLocal(INSERT_ALL_VALUES) to fill it up

  3) DMClone(dm) and give it a Section with no constraints

  4) Use DMLocalToGlobal() into the new DM

  5) VecView() that new global vector

  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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150423/8c7c99e5/attachment.html>