Hi, I am trying to write a second order upwind finite volume scheme on unstructured grid using DMPlex. I am performing gradient reconstruction using the "DMPlexReconstructGradientsFVM" function as well as a "ComputeGradient" function that I have written myself, in which, I loop over faces and add the gradient contribution to neighboring cells of the face (similar to what is done in "DMPlexReconstructGradients_Internal") . The results don't seem correct when using the "DMPlexReconstructGradientsFVM" function where as when using my own function, the second order is achieved. To investigate, I looked at the gradients computed by these two functions, and noticed that "DMPlexReconstructGradientsFVM" assigns a zero value to the gradient of cells with less than three neighbors, for example, boundary cells (which I assume is because leastsquares requires over-constrained system of equations and hence, atleast three neighbors), where as, I am not considering this factor in my "ComputeGradient" function and adding the contribution to neighboring cells for all interior and partition faces. Other than this, both the functions seem to compute similar values for gradients on all the other cells and I suspect, zero values for certain cells might be causing the problem. I wanted to know if I am missing out on some crucial step in the code which would enable the "DMPlexReconstructGradientsFVM" function to compute gradients for cells with lesser number of neighbors, and also if there is a way to augment the leastsquares stencil, i.e. to use, say, vertex neighbors instead of face neighbors to reconstruct gradients. Also, kindly let me know if my analysis of the problem is wrong and if there might be some other issue.
I am attaching my code for your reference. You can find the use of both
"DMPlexReconstructGradientsFVM" function and "ComputeGradient" function
(just uncomment the one you want to use) inside "ComputeResidual" function.
I am also attaching the Gmsh file for the mesh I am using. You can run the
code as follows:
mpiexec -n <NumProc> ./convect -mesh square_tri.msh
Please let me know if you need anything else.
Thanks and Regards,
Shashwat
static char help[] = "Linear advection using first order FVM\n";
#include <petsc.h>
PETSC_STATIC_INLINE
void Waxpy(PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
{
PetscInt d; for (d = 0; d < dim; ++d) w[d] = a*x[d] + y[d];
}
PETSC_STATIC_INLINE
PetscScalar Dot(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
{
PetscScalar sum = 0.0; PetscInt d; for (d = 0; d < dim; ++d) sum += x[d]*y[d]; return sum;
}
PETSC_STATIC_INLINE
PetscReal Norm(PetscInt dim, const PetscScalar *x)
{
return PetscSqrtReal(PetscAbsScalar(Dot(dim,x,x)));
}
typedef struct
{
PetscScalar *X0, *V0; // cell data
PetscScalar *X1, *V1, *N1; // face data
PetscScalar *Gl, *Gr; // gradient contributions
PetscInt ncell, nface;
DM dmGrad;
} MeshInfo;
typedef struct
{
MeshInfo info;
PetscInt vtkInterval;
// index sets for boundary, interior and partition faces
IS is_bface, is_iface, is_pface;
} AppCtx;
// advection speed
void advection_speed(const PetscReal *X, PetscReal *speed)
{
PetscReal x = X[0];
PetscReal y = X[1];
speed[0] = -y;
speed[1] = x;
}
// initial condition
PetscReal initial_condition(const PetscReal *X)
{
PetscReal x = X[0];
PetscReal y = X[1];
PetscReal r2 = PetscPowReal(x-0.5,2) + PetscPowReal(y,2);
return PetscExpReal(-50.0*r2);
}
// upwind numerical flux
PetscReal numerical_flux(const PetscReal a, const PetscReal ul, const PetscReal ur)
{
return (a > 0.0) ? a*(ul) : a*(ur);
}
// construct label to differentiate real and ghost cells
// LabelValue = 1 => Real Cells
// 0 => Ghost Cells (both boundary and partition ghosts)
PetscErrorCode ConstructCellLabels(DM dm, DMLabel *rCell)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscInt c, cStart, cEnd; // cells
PetscInt gStart, gEnd; // ghost cells (boundary ghosts)
// get ranges for all cells and ghost (boundary) cells
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
ierr = DMPlexGetGhostCellStratum(dm, &gStart, &gEnd); CHKERRQ(ierr);
// loop over real cells (includes partition ghosts)
for(c=cStart; c<gStart; ++c)
{
PetscInt lVal; // for partition ghosts
ierr = DMGetLabelValue(dm, "ghost", c, &lVal); CHKERRQ(ierr);
if(lVal != 2) // not a partition ghost
{
ierr = DMLabelSetValue(*rCell, c, 1); CHKERRQ(ierr);
}
else
{
ierr = DMLabelSetValue(*rCell, c, 0); CHKERRQ(ierr);
}
}
// mark boundary ghost cells
for(c=gStart; c<gEnd; ++c)
{
ierr = DMLabelSetValue(*rCell, c, 0); CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
// construct label for internal and boundary faces with index sets
// LabelValue = 0 => partition boundary face
// = 1 => real boundary face
// = 2 => real interior face
// = -1 => ghost face
PetscErrorCode ConstructFaceLabels(DM dm, DMLabel *faces)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscInt f, fStart, fEnd; // faces
PetscInt ss;
const PetscInt *supp; // support
ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd); CHKERRQ(ierr);
// loop over faces and assign label
for(f=fStart; f<fEnd; ++f)
{
ierr = DMPlexGetSupportSize(dm, f, &ss); CHKERRQ(ierr);
ierr = DMPlexGetSupport(dm, f, &supp); CHKERRQ(ierr);
if(ss == 1)
{
PetscInt label;
ierr = DMGetLabelValue(dm, "cells", supp[0], &label); CHKERRQ(ierr);
// check if the face has real neighboring cell
if(label == 1) // real boundary face
{
ierr = DMLabelSetValue(*faces, f, 1); CHKERRQ(ierr);
}
}
else if(ss == 2)
{
PetscInt label1, label2;
ierr = DMGetLabelValue(dm, "cells", supp[0], &label1); CHKERRQ(ierr);
ierr = DMGetLabelValue(dm, "cells", supp[1], &label2); CHKERRQ(ierr);
if(label1 == 1 && label2 == 1) // both real cells
{
ierr = DMLabelSetValue(*faces, f, 2); CHKERRQ(ierr);
}
else if(!(label1 != 1 && label2 != 1)) // atleast one real cell
{
PetscInt glabel; // ghost label (to check if the neighbor is partition ghost)
ierr = DMGetLabelValue(dm, "ghost", supp[1], &glabel); CHKERRQ(ierr);
if(glabel == 2) // partition ghost
{
ierr = DMLabelSetValue(*faces, f, 0); CHKERRQ(ierr); // partition boundary
}
else
{
ierr = DMLabelSetValue(*faces, f, 1); CHKERRQ(ierr); // real boundary
}
}
}
}
PetscFunctionReturn(0);
}
// set initial condition for the given dm and vector
PetscErrorCode SetIC(DM dm, Vec U)
{
PetscErrorCode ierr;
PetscScalar *u;
PetscFunctionBegin;
PetscInt c, cStart, cEnd; // cells
PetscReal area, centroid[3], normal[3]; // geometric data
// get cell stratum owned by processor
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
// get array for U
ierr = VecGetArray(U, &u);
// loop over cells and assign values
for(c=cStart; c<cEnd; ++c)
{
PetscInt label;
ierr = DMGetLabelValue(dm, "cells", c, &label); CHKERRQ(ierr);
// write into Global vector if the cell is a real cell
if(label == 1)
{
PetscReal X[2]; // cell centroid
ierr = DMPlexComputeCellGeometryFVM(dm, c, &area, centroid, normal); CHKERRQ(ierr);
X[0] = centroid[0]; X[1] = centroid[1];
u[c] = initial_condition(X);
}
}
ierr = VecRestoreArray(U, &u);
PetscFunctionReturn(0);
}
static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
{
PetscFunctionBegin;
PetscErrorCode ierr;
ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(*viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(*viewer, PETSC_VIEWER_VTK_VTU); CHKERRQ(ierr);
ierr = PetscViewerFileSetName(*viewer, filename);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
// compute gradients on each cell
PetscErrorCode ComputeGradient(DM dm, Vec U, Vec grad, MeshInfo *info)
{
PetscErrorCode ierr;
PetscFunctionBegin;
PetscInt f, fStart, fEnd; // faces
PetscInt cStart, cEnd; // cells
const PetscInt *supp, *face;
ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd); CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
IS is;
PetscInt fSize;
Vec Ul;
ierr = DMGetLocalVector(dm, &Ul); CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Ul); CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Ul); CHKERRQ(ierr);
ierr = VecSet(grad, 0.0); CHKERRQ(ierr);
PetscScalar *UArr, *gArr;
ierr = VecGetArray(Ul, &UArr); CHKERRQ(ierr);
ierr = VecGetArray(grad, &gArr); CHKERRQ(ierr);
// interior face
ierr = DMGetStratumIS(dm, "faces", 2, &is); CHKERRQ(ierr);
ierr = ISGetLocalSize(is, &fSize); CHKERRQ(ierr);
ierr = ISGetIndices(is, &face); CHKERRQ(ierr);
for(f=0; f<fSize; ++f)
{
ierr = DMPlexGetSupport(dm, face[f], &supp); CHKERRQ(ierr);
PetscReal ul = UArr[supp[0]], ur = UArr[supp[1]];
// left cell
gArr[2*supp[0]] += info->Gl[2*(face[f]-fStart)]*(ur - ul);
gArr[2*supp[0]+1] += info->Gl[2*(face[f]-fStart)+1]*(ur - ul);
// right cell
gArr[2*supp[1]] -= info->Gr[2*(face[f]-fStart)]*(ur - ul);
gArr[2*supp[1]+1] -= info->Gr[2*(face[f]-fStart)+1]*(ur - ul);
}
// partition face
// interior face
ierr = DMGetStratumIS(dm, "faces", 0, &is); CHKERRQ(ierr);
if(is)
{
ierr = ISGetLocalSize(is, &fSize); CHKERRQ(ierr);
ierr = ISGetIndices(is, &face); CHKERRQ(ierr);
for(f=0; f<fSize; ++f)
{
ierr = DMPlexGetSupport(dm, face[f], &supp); CHKERRQ(ierr);
PetscReal ul = UArr[supp[0]], ur = UArr[supp[1]];
// left cell
gArr[2*supp[0]] += info->Gl[2*(face[f]-fStart)]*(ur - ul);
gArr[2*supp[0]+1] += info->Gl[2*(face[f]-fStart)+1]*(ur - ul);
}
}
ierr = VecRestoreArray(Ul, &UArr); CHKERRQ(ierr);
ierr = VecRestoreArray(grad, &gArr); CHKERRQ(ierr);
ierr = DMRestoreLocalVector(dm, &Ul); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
// function to compute residual
PetscErrorCode ComputeResidual(TS ts, PetscReal time, Vec U, Vec R, void *ctx)
{
PetscErrorCode ierr;
AppCtx *user = (AppCtx *) ctx;
MeshInfo *info = &user->info;
PetscFunctionBegin;
DM dm;
Vec grad;
Vec Ul, locGrad; // local vectors
// get dm from ts
ierr = TSGetDM(ts, &dm); CHKERRQ(ierr);
// get local vector from dm
ierr = DMGetLocalVector(dm, &Ul); CHKERRQ(ierr);
// initialize R
ierr = VecSet(R, 0.0); CHKERRQ(ierr);
// copy values from Global to Local vector (including ghost values)
ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, Ul); CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, Ul); CHKERRQ(ierr);
// compute grandients using least squares
ierr = DMCreateGlobalVector(info->dmGrad, &grad); CHKERRQ(ierr);
ierr = DMPlexReconstructGradientsFVM(dm, Ul, grad); CHKERRQ(ierr);
//ierr = ComputeGradient(dm, U, grad, info); CHKERRQ(ierr);
ierr = DMCreateLocalVector(info->dmGrad, &locGrad); CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(info->dmGrad, grad, INSERT_VALUES, locGrad); CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(info->dmGrad, grad, INSERT_VALUES, locGrad); CHKERRQ(ierr);
// get array for local gradients
PetscReal *gradArr, *gradComp;
ierr = VecGetArray(locGrad, &gradArr); CHKERRQ(ierr);
// get solution and residual arrays
PetscReal *Rarr, *Uarr;
ierr = VecGetArray(Ul, &Uarr); CHKERRQ(ierr);
ierr = VecGetArray(R, &Rarr); CHKERRQ(ierr);
PetscInt f, fStart, fEnd, fSize; // faces
PetscInt cStart, cEnd; // cells
const PetscInt *supp, *face;
PetscReal ul, ur, flux, norm_speed, speed[2], Xfc[2];
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd); CHKERRQ(ierr);
// loop over faces
// real boundary faces
ierr = ISGetLocalSize(user->is_bface, &fSize); CHKERRQ(ierr);
ierr = ISGetIndices(user->is_bface, &face); CHKERRQ(ierr); // get indices
for(f=0; f<fSize; ++f)
{
// get support
ierr = DMPlexGetSupport(dm, face[f], &supp); CHKERRQ(ierr);
// calculate speed
advection_speed(&info->X1[2*(face[f]-fStart)], speed);
norm_speed = Dot(2, speed, &info->N1[2*(face[f]-fStart)]);
// get gradient for cell
ierr = DMPlexPointLocalRead(info->dmGrad, supp[0], gradArr, &gradComp); CHKERRQ(ierr);
// get vector from cell centroid to face centroid
Xfc[0] = info->X1[2*(face[f]-fStart)] - info->X0[2*(supp[0]-cStart)];
Xfc[1] = info->X1[2*(face[f]-fStart)+1] - info->X0[2*(supp[0]-cStart)+1];
ul = Uarr[supp[0]] + Dot(2, gradComp, Xfc);
ur = 0; // Dirichlet BC
flux = numerical_flux(norm_speed, ul, ur);
Rarr[supp[0]] -= flux * info->V1[face[f]-fStart] / info->V0[supp[0]-cStart];
}
ierr = ISRestoreIndices(user->is_bface, &face); CHKERRQ(ierr); // restore index set
if(user->is_pface) // if there are partition faces
{
ierr = ISGetLocalSize(user->is_pface, &fSize); CHKERRQ(ierr);
ierr = ISGetIndices(user->is_pface, &face); CHKERRQ(ierr); // get indices
for(f=0; f<fSize; ++f)
{
// get support
ierr = DMPlexGetSupport(dm, face[f], &supp); CHKERRQ(ierr);
// calculate speed
advection_speed(&info->X1[2*(face[f]-fStart)], speed);
norm_speed = Dot(2, speed, &info->N1[2*(face[f]-fStart)]);
// get gradient for cell
ierr = DMPlexPointLocalRead(info->dmGrad, supp[0], gradArr, &gradComp); CHKERRQ(ierr);
// get vector from cell centroid to face centroid
Xfc[0] = info->X1[2*(face[f]-fStart)] - info->X0[2*(supp[0]-cStart)];
Xfc[1] = info->X1[2*(face[f]-fStart)+1] - info->X0[2*(supp[0]-cStart)+1];
ul = Uarr[supp[0]] + Dot(2, gradComp, Xfc);
// get gradient for cell
ierr = DMPlexPointLocalRead(info->dmGrad, supp[1], gradArr, &gradComp); CHKERRQ(ierr);
// get vector from cell centroid to face centroid
Xfc[0] = info->X1[2*(face[f]-fStart)] - info->X0[2*(supp[1]-cStart)];
Xfc[1] = info->X1[2*(face[f]-fStart)+1] - info->X0[2*(supp[1]-cStart)+1];
ur = Uarr[supp[1]] + Dot(2, gradComp, Xfc);
flux = numerical_flux(norm_speed, ul, ur);
Rarr[supp[0]] -= flux * info->V1[face[f]-fStart] / info->V0[supp[0]-cStart];
}
ierr = ISRestoreIndices(user->is_pface, &face); CHKERRQ(ierr);
}
// interior faces
ierr = ISGetLocalSize(user->is_iface, &fSize); CHKERRQ(ierr);
ierr = ISGetIndices(user->is_iface, &face); CHKERRQ(ierr); // get indices
for(f=0; f<fSize; ++f)
{
// get support
ierr = DMPlexGetSupport(dm, face[f], &supp); CHKERRQ(ierr);
advection_speed(&info->X1[2*(face[f]-fStart)], speed);
norm_speed = Dot(2, speed, &info->N1[2*(face[f]-fStart)]);
// get gradient for cell
ierr = DMPlexPointLocalRead(info->dmGrad, supp[0], gradArr, &gradComp); CHKERRQ(ierr);
// get vector from cell centroid to face centroid
Xfc[0] = info->X1[2*(face[f]-fStart)] - info->X0[2*(supp[0]-cStart)];
Xfc[1] = info->X1[2*(face[f]-fStart)+1] - info->X0[2*(supp[0]-cStart)+1];
ul = Uarr[supp[0]] + Dot(2, gradComp, Xfc);
// get gradient for cell
ierr = DMPlexPointLocalRead(info->dmGrad, supp[1], gradArr, &gradComp); CHKERRQ(ierr);
// get vector from cell centroid to face centroid
Xfc[0] = info->X1[2*(face[f]-fStart)] - info->X0[2*(supp[1]-cStart)];
Xfc[1] = info->X1[2*(face[f]-fStart)+1] - info->X0[2*(supp[1]-cStart)+1];
ur = Uarr[supp[1]] + Dot(2, gradComp, Xfc);
flux = numerical_flux(norm_speed, ul, ur);
Rarr[supp[0]] -= flux * info->V1[face[f]-fStart] / info->V0[supp[0]-cStart];
Rarr[supp[1]] += flux * info->V1[face[f]-fStart] / info->V0[supp[1]-cStart];
}
ierr = ISRestoreIndices(user->is_iface, &face); CHKERRQ(ierr); // restore index set
// cleanup
ierr = VecRestoreArray(locGrad, &gradArr); CHKERRQ(ierr);
ierr = VecRestoreArray(Ul, &Uarr); CHKERRQ(ierr);
ierr = VecRestoreArray(R, &Rarr); CHKERRQ(ierr);
ierr = DMRestoreLocalVector(dm, &Ul); CHKERRQ(ierr);
ierr = DMRestoreLocalVector(info->dmGrad, &locGrad); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(info->dmGrad, &grad); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
// monitor function for the TS context
PetscErrorCode MonitorTS(TS ts, PetscInt step, PetscReal time, Vec U, void *ctx)
{
PetscErrorCode ierr;
static PetscInt counter = 0;
AppCtx *user = (AppCtx *) ctx;
PetscInt vtkInterval = user->vtkInterval;
char filename[PETSC_MAX_PATH_LEN];
PetscViewer viewer;
PetscReal tf;
DM dm;
PetscFunctionBegin;
ierr = PetscPrintf(PETSC_COMM_WORLD, "Iter: %d, t = %f\n", step, time); CHKERRQ(ierr);
ierr = TSGetMaxTime(ts, &tf); CHKERRQ(ierr);
if(step % vtkInterval == 0 || time == tf)
{
ierr = TSGetDM(ts, &dm); CHKERRQ(ierr);
ierr = PetscSNPrintf(filename, sizeof(filename),
"sol-%03D.vtu", counter); CHKERRQ(ierr);
ierr = OutputVTK(dm, filename, &viewer); CHKERRQ(ierr);
ierr = VecView(U, viewer); CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Wrote %s\n", filename); CHKERRQ(ierr);
++counter;
}
PetscFunctionReturn(0);
}
// setup Finite Volume object
PetscErrorCode SetUpFV(PetscFV fv)
{
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscFVSetType(fv, PETSCFVLEASTSQUARES); CHKERRQ(ierr);
ierr = PetscFVSetNumComponents(fv, 1); CHKERRQ(ierr);
ierr = PetscFVSetSpatialDimension(fv, 2); CHKERRQ(ierr);
ierr = PetscFVSetFromOptions(fv); CHKERRQ(ierr);
ierr = PetscFVSetUp(fv); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
// create geometric data for cells and faces
PetscErrorCode CreateMeshInfo(DM dm, MeshInfo *info)
{
PetscErrorCode ierr;
PetscFunctionBegin;
Vec cellGeom, faceGeom;
DM dmFace, dmCell;
PetscObject fv;
ierr = DMGetField(dm, 0, NULL, &fv); CHKERRQ(ierr);
const PetscScalar *fgeom, *cgeom; // arrays corresponding to vectors
PetscFVCellGeom *cg; // struct with cell geometry information
PetscFVFaceGeom *fg; // struct with face geometry information
ierr = DMPlexGetDataFVM(dm, (PetscFV) fv, &cellGeom, &faceGeom, &info->dmGrad); CHKERRQ(ierr);
ierr = VecGetDM(faceGeom, &dmFace); CHKERRQ(ierr);
ierr = VecGetDM(cellGeom, &dmCell); CHKERRQ(ierr);
ierr = VecGetArrayRead(faceGeom, &fgeom); CHKERRQ(ierr);
ierr = VecGetArrayRead(cellGeom, &cgeom); CHKERRQ(ierr);
PetscInt p, pStart, pEnd; // points on DM
// cell info
ierr = DMPlexGetHeightStratum(dm, 0, &pStart, &pEnd); CHKERRQ(ierr);
info->ncell = pEnd - pStart;
ierr = PetscMalloc2(2*info->ncell, &info->X0, info->ncell, &info->V0); CHKERRQ(ierr);
for(p=0; p<pEnd-pStart; ++p)
{
ierr = DMPlexPointLocalRead(dmCell, p + pStart, cgeom, &cg); CHKERRQ(ierr);
info->X0[2*p] = cg->centroid[0];
info->X0[2*p+1] = cg->centroid[1];
info->V0[p] = cg->volume;
}
// face info
ierr = DMPlexGetHeightStratum(dm, 1, &pStart, &pEnd); CHKERRQ(ierr);
info->nface = pEnd - pStart;
ierr = PetscMalloc3(2*info->nface, &info->X1, 2*info->nface, &info->N1,
info->nface, &info->V1); CHKERRQ(ierr);
ierr = PetscMalloc2(2*info->nface, &info->Gl, 2*info->nface, &info->Gr); CHKERRQ(ierr);
for(p=0; p<pEnd-pStart; ++p)
{
ierr = DMPlexPointLocalRead(dmFace, p + pStart, fgeom, &fg); CHKERRQ(ierr);
info->X1[2*p] = fg->centroid[0];
info->X1[2*p+1] = fg->centroid[1];
info->V1[p] = Norm(2, fg->normal);
info->N1[2*p] = fg->normal[0]/info->V1[p]; // normal is scaled
info->N1[2*p+1] = fg->normal[1]/info->V1[p];
info->Gl[2*p] = fg->grad[0][0];
info->Gl[2*p+1] = fg->grad[0][1];
info->Gr[2*p] = fg->grad[1][0];
info->Gr[2*p+1] = fg->grad[1][1];
}
// cleanup
ierr = VecRestoreArrayRead(faceGeom, &fgeom); CHKERRQ(ierr);
ierr = VecRestoreArrayRead(cellGeom, &cgeom); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
// destroy mesh info
PetscErrorCode DestroyMeshInfo(MeshInfo *info)
{
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscFree2(info->X0, info->V0); CHKERRQ(ierr);
ierr = PetscFree3(info->X1, info->N1, info->V1); CHKERRQ(ierr);
ierr = PetscFree2(info->Gl, info->Gr); CHKERRQ(ierr);
PetscFunctionReturn(0);
}
int main(int argc, char **argv)
{
PetscErrorCode ierr;
DM dm, dmDist = NULL;
TS ts;
PetscSection sec;
PetscInt dim = 2, numFields = 1, overlap = 1, numBC, i;
PetscInt numComp[1]; // number of components per field
PetscInt numDof[3];
PetscInt bcField[1]; // number of BCs for each field
IS bcPointsIS[1]; // contains boundary faces
PetscBool interpolate = PETSC_TRUE, status;
Vec U, R; // solution and residual
PetscMPIInt rank, size;
PetscReal lower[3], upper[3]; // lower left and upper right coordinates of domain
PetscInt cells[2];
char filename[PETSC_MAX_PATH_LEN];
AppCtx user;
// define domain properties
cells[0] = 100; cells[1] = 100;
lower[0] = -1; lower[1] = -1; lower[2] = 0;
upper[0] = 1; upper[1] = 1; upper[2] = 0;
ierr = PetscInitialize(&argc, &argv, (char*)0, help); CHKERRQ(ierr);
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);
// get number of cells in x and y direction
ierr = PetscOptionsGetInt(NULL, NULL, "-nx", &cells[0], &status); CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL, NULL, "-ny", &cells[1], &status); CHKERRQ(ierr);
// set time parameters
PetscReal tf = 2.0*M_PI, dt = 0.001;
user.vtkInterval = 100;
ierr = PetscOptionsGetReal(NULL, NULL, "-tf", &tf, &status); CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL, NULL, "-vtk", &user.vtkInterval, &status); CHKERRQ(ierr);
// get mesh filename
ierr = PetscOptionsGetString(NULL, NULL, "-mesh", filename, PETSC_MAX_PATH_LEN, &status);
if(status) // gmsh file provided by user
{
char file[PETSC_MAX_PATH_LEN];
ierr = PetscStrcpy(file, filename); CHKERRQ(ierr);
ierr = PetscSNPrintf(filename, sizeof filename,"./%s", file); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Reading gmsh %s ... ", file); CHKERRQ(ierr);
ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, filename, interpolate, &dm); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); CHKERRQ(ierr);
}
else
{
// create the mesh
ierr = PetscPrintf(PETSC_COMM_WORLD, "Generating mesh ... "); CHKERRQ(ierr);
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE,
cells, lower, upper, NULL, interpolate, &dm); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Done\n"); CHKERRQ(ierr);
}
// distribute mesh over processes;
ierr = DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE); CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, NULL, &dmDist); CHKERRQ(ierr);
if(dmDist)
{
ierr = DMDestroy(&dm); CHKERRQ(ierr);
dm = dmDist;
}
// print mesh information
ierr = PetscPrintf(PETSC_COMM_WORLD, "overlap: %d, "
"distributed among %d processors\n",
overlap, size); CHKERRQ(ierr);
// set up FVM
PetscFV fv;
ierr = PetscFVCreate(PETSC_COMM_WORLD, &fv); CHKERRQ(ierr);
ierr = SetUpFV(fv); CHKERRQ(ierr);
ierr = DMSetNumFields(dm, numFields); CHKERRQ(ierr);
ierr = DMSetField(dm, 0, NULL, (PetscObject) fv); CHKERRQ(ierr);
ierr = DMCreateDS(dm); CHKERRQ(ierr);
ierr = PetscFVDestroy(&fv); CHKERRQ(ierr);
// construct ghost cells
PetscInt nGhost; // number of ghost cells
DM dmG; // DM with ghost cells
ierr = DMPlexConstructGhostCells(dm, NULL, &nGhost, &dmG); CHKERRQ(ierr);
if(dmG)
{
ierr = DMDestroy(&dm); CHKERRQ(ierr);
dm = dmG;
}
// construct label for real cells
DMLabel rCell;
const char *cLabel = "cells";
ierr = DMCreateLabel(dm, cLabel); CHKERRQ(ierr);
ierr = DMGetLabel(dm, cLabel, &rCell); CHKERRQ(ierr);
ierr = ConstructCellLabels(dm, &rCell); CHKERRQ(ierr);
// construct face labels
DMLabel faces;
const char *fLabel = "faces";
ierr = DMCreateLabel(dm, fLabel); CHKERRQ(ierr);
ierr = DMGetLabel(dm, fLabel, &faces); CHKERRQ(ierr);
ierr = ConstructFaceLabels(dm, &faces); CHKERRQ(ierr);
// create mesh info
ierr = CreateMeshInfo(dm, &user.info); CHKERRQ(ierr);
// get index sets for faces
ierr = DMGetStratumIS(dm, "faces", 0, &user.is_pface); CHKERRQ(ierr);
ierr = DMGetStratumIS(dm, "faces", 1, &user.is_bface); CHKERRQ(ierr);
ierr = DMGetStratumIS(dm, "faces", 2, &user.is_iface); CHKERRQ(ierr);
// create scalar field for solution u
numComp[0] = 1;
for(i=0; i<numFields*(dim+1); ++i) numDof[i] = 0;
// define u at cell centers
numDof[0*(dim+1)+dim] = 1;
// set BC
numBC = 1;
bcField[0] = 0; // Dirichlet BC on u at boundary
ierr = DMGetStratumIS(dm, "faces", 1, &bcPointsIS[0]); CHKERRQ(ierr);
// create a PetscSection with this layout
ierr = DMPlexCreateSection(dm, NULL, numComp, numDof, numBC,
bcField, NULL, bcPointsIS, NULL, &sec); CHKERRQ(ierr);
ierr = ISDestroy(&bcPointsIS[0]); CHKERRQ(ierr);
ierr = PetscSectionSetFieldName(sec, 0, "u"); CHKERRQ(ierr);
// set the section for dm
ierr = DMSetLocalSection(dm, sec); CHKERRQ(ierr);
ierr = DMSetUp(dm); CHKERRQ(ierr);
ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);
// get global vector for solution and set IC
ierr = DMCreateGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) U, "u"); CHKERRQ(ierr);
ierr = SetIC(dm, U); CHKERRQ(ierr);
ierr = VecDuplicate(U, &R); CHKERRQ(ierr);
// setup the TS context for time stepping and solve
ierr = TSCreate(PETSC_COMM_WORLD, &ts); CHKERRQ(ierr);
ierr = TSSetProblemType(ts, TS_NONLINEAR); CHKERRQ(ierr);
ierr = TSSetRHSFunction(ts, R, ComputeResidual, &user); CHKERRQ(ierr);
ierr = TSSetMaxTime(ts, tf); CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP); CHKERRQ(ierr);
ierr = TSMonitorSet(ts, MonitorTS, &user, NULL); CHKERRQ(ierr);
ierr = TSSetDM(ts, dm); CHKERRQ(ierr);
ierr = TSSetTimeStep(ts, dt); CHKERRQ(ierr);
ierr = TSRKSetType(ts, TSRK2A); CHKERRQ(ierr);
ierr = TSSetType(ts, TSRK); CHKERRQ(ierr);
ierr = TSSetSolution(ts, U); CHKERRQ(ierr);
ierr = TSSetFromOptions(ts); CHKERRQ(ierr);
ierr = TSSolve(ts, U); CHKERRQ(ierr);
// cleanup
ierr = DestroyMeshInfo(&user.info); CHKERRQ(ierr);
ierr = ISDestroy(&user.is_pface); CHKERRQ(ierr);
ierr = ISDestroy(&user.is_bface); CHKERRQ(ierr);
ierr = ISDestroy(&user.is_iface); CHKERRQ(ierr);
ierr = VecDestroy(&R); CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm, &U); CHKERRQ(ierr);
ierr = PetscSectionDestroy(&sec); CHKERRQ(ierr);
ierr = TSDestroy(&ts); CHKERRQ(ierr);
ierr = DMDestroy(&dm); CHKERRQ(ierr);
ierr = PetscFinalize(); CHKERRQ(ierr);
return ierr;
}
square_tri.geo
Description: Binary data
