Skip to content

Commit

Permalink
Solve projection problems with two elements. [#1 status:resolved]
Browse files Browse the repository at this point in the history 
Two main weaknesses:

* Uses -snes_mf for projection, we should also assemble a matrix and be
  able to apply the Jacobian matrix-free.

* Discretization error is not yet checked, but evaluation at quadrature
  points demonstrates spectral convergence so the projection should as
  well.

Signed-off-by: Jed Brown <jed@59A2.org>
  • Loading branch information
@jedbrown
jedbrown committed Nov 21, 2008
1 parent 87ffc4e commit d022dba
Show file tree
Hide file tree
Showing 16 changed files with 467 additions and 89 deletions.
4 changes: 4 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,10 @@ if (Dohp_USE_DEBUG)
add_definitions (-DdUSE_DEBUG=1 -g3)
endif (Dohp_USE_DEBUG)

if (Dohp_BUILD_TESTS)
enable_testing ()
endif (Dohp_BUILD_TESTS)

if (PEDANTIC_WARNINGS)
set (DEFAULT_PEDANTIC_FLAGS "-pedantic -Wall -Wextra -Wundef -Wshadow -Wpointer-arith -Wbad-function-cast -Wcast-align -Wwrite-strings -Wconversion -Wlogical-op -Wsign-compare -Waggregate-return -Wstrict-prototypes -Wmissing-prototypes -Wmissing-declarations -Wredundant-decls -Wnested-externs -Winline -Wno-long-long -Wmissing-format-attribute -Wmissing-noreturn -Wpacked -Wdisabled-optimization -Wmultichar -Wformat-nonliteral -Wformat-security -Wformat-y2k -Wendif-labels -Wdeclaration-after-statement -Wold-style-definition -Winvalid-pch -Wmissing-field-initializers -Wvariadic-macros -Wunsafe-loop-optimizations -Wvolatile-register-var -Wstrict-aliasing -funit-at-a-time -Wno-sign-conversion")
#set (DEFAULT_PEDANTIC_FLAGS "-Wunreachable-code -Wfloat-equal -Wc++-compat")
Expand Down
6 changes: 6 additions & 0 deletions include/dohpfs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,12 @@ EXTERN dErr dFSCreateGlobalVector(dFS,Vec*);
EXTERN dErr dFSGlobalToExpandedBegin(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSGlobalToExpandedEnd(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToGlobal(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToGlobalBegin(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSExpandedToGlobalEnd(dFS,Vec,InsertMode,Vec);
EXTERN dErr dFSGetElements(dFS,dInt*,dInt**,s_dRule**,s_dEFS**,dInt**,dReal(**)[3]);
EXTERN dErr dFSRestoreElements(dFS,dInt*,dInt**,s_dRule**,s_dEFS**,dInt**,dReal(**)[3]);
EXTERN dErr dFSGetWorkspace(dFS,dReal(**)[3],dReal(**)[3][3],dReal**,dScalar**,dScalar**,dScalar**,dScalar**);
EXTERN dErr dFSRestoreWorkspace(dFS,dReal(**)[3],dReal(**)[3][3],dReal**,dScalar**,dScalar**,dScalar**,dScalar**);
EXTERN dErr dFSMatSetValuesExpanded(dFS,Mat,dInt,const dInt[],dInt,const dInt[],const dScalar[],InsertMode);
EXTERN dErr dFSGetMatrix(dFS,const MatType,Mat*);
EXTERN dErr dFSBuildSpace(dFS);
Expand Down
3 changes: 2 additions & 1 deletion include/dohpjacobi.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,7 +151,8 @@ EXTERN dErr dJacobiGetEFS(dJacobi,dInt,const dEntTopology[],const dInt[],dRule,d
EXTERN dErr dRuleView(dRule rule,PetscViewer);
EXTERN dErr dRuleGetSize(dRule rule,dInt *dim,dInt *nnodes);
EXTERN dErr dRuleGetNodeWeight(dRule rule,dReal *coord,dReal *weight);
EXTERN dErr dRuleGetTensorNodeWeight(dRule rule,dInt *dim,dInt *nnodes,const dReal **coord,const dReal **weight);
EXTERN dErr dRuleGetTensorNodeWeight(dRule rule,dInt *dim,dInt *nnodes,const dReal *coord[],const dReal *weight[]);
EXTERN dErr dRuleComputeGeometry(dRule rule,const dReal vtx[restrict][3],dReal[restrict][3],dReal jinv[restrict][3][3],dReal jdet[restrict]);

EXTERN dErr dEFSView(dEFS efs,PetscViewer viewer);
EXTERN dErr dEFSGetSizes(dEFS efs,dInt*,dInt *inodes,dInt *total);
Expand Down
2 changes: 2 additions & 0 deletions include/dohpmesh.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,8 @@ EXTERN dErr dMeshGetStatus(dMesh,dInt,const dMeshEH[],dEntStatus[]);
EXTERN dErr dMeshGetTopo(dMesh,dInt,const dMeshEH[],dEntTopology[]);
EXTERN dErr dMeshGetAdjacency(dMesh,dMeshESH,struct dMeshAdjacency*);
EXTERN dErr dMeshRestoreAdjacency(dMesh,dMeshESH,struct dMeshAdjacency*);
EXTERN dErr dMeshGetVertexCoords(dMesh,dInt,const dMeshEH[],dInt**,dReal(**)[3]);
EXTERN dErr dMeshRestoreVertexCoords(dMesh,dInt,const dMeshEH[],dInt**,dReal(**)[3]);

PETSC_EXTERN_CXX_END
#endif
1 change: 1 addition & 0 deletions include/dohptype.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -117,6 +117,7 @@ typedef unsigned char dEntStatus;
#define dValidRealPointer(p,a) dValidPointerNamedSpecific((p),dReal,"dReal",(a))

#define dStrlen(s,l) PetscStrlen((s),(l))
static inline dErr dObjectGetComm(dObject obj,MPI_Comm *comm) { return PetscObjectGetComm(obj,comm); }

static inline dInt dMaxInt(dInt a,dInt b) { return (a > b) ? a : b; }
static inline dInt dMinInt(dInt a,dInt b) { return (a < b) ? a : b; }
Expand Down
55 changes: 34 additions & 21 deletions include/private/fsimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,14 @@ struct _p_dFSBoundary {
struct _p_dFSBoundary *next;
};

typedef struct {
dInt Q;
dReal (*q)[3];
dReal (*jinv)[3][3];
dReal *jw;
dScalar *u,*v,*du,*dv;
} s_dFSWorkspace;

struct _dFSOps {
DMOPS(dFS)
dErr (*setfromoptions)(dFS);
Expand All @@ -28,27 +36,32 @@ struct _dFSOps {
struct _p_dFS {
PETSCHEADER(struct _dFSOps);
DMHEADER
dMesh mesh;
dMeshTag degreetag,ruletag; /**< tags on regions */
dMeshESH active; /**< regions that will be part of this space */
dFSBoundary bdylist;
dQuotient quotient;
dJacobi jacobi;
dTruth spacebuilt;
dTruth assemblefull; /**< Use full order constraints for assembly */
Sliced sliced;
dInt nbdofs;
dInt n,N; /**< length of the owned and global vectors */
dInt nlocal; /**< number of owned+ghost dofs on this process */
dInt rstart; /**< global offset of first owned dof */
dInt m; /**< Number of expanded dofs */
dInt D; /**< Number of dofs per (non-boundary) node */
s_dRule *rule; /**< Integration rule */
s_dEFS *efs; /**< Element function space, defined for all entities */
Mat C; /**< full-order constraint matrix (element dofs to local numbering) */
Mat Cp; /**< preconditioning constraint matrix (element dofs to local numbering, as sparse as possible) */
Vec weight; /**< Vector in global space, used to compensate for overcounting after local to global */
void *data;
dMesh mesh;
dMeshTag degreetag,ruletag; /**< tags on regions */
dMeshESH active; /**< regions that will be part of this space */
dFSBoundary bdylist;
dQuotient quotient;
dJacobi jacobi;
dTruth spacebuilt;
dTruth assemblefull; /**< Use full order constraints for assembly */
Sliced sliced;
dInt nbdofs;
dInt n,N; /**< length of the owned and global vectors */
dInt nlocal; /**< number of owned+ghost dofs on this process */
dInt rstart; /**< global offset of first owned dof */
dInt m; /**< Number of expanded dofs */
dInt D; /**< Number of dofs per (non-boundary) node */
dInt nelem;
dInt *off; /**< Offset of element dofs in expanded vector */
s_dRule *rule; /**< Integration rule */
s_dEFS *efs; /**< Element function space, defined for all entities */
dInt *vtxoff;
dReal (*vtx)[3];
Mat C; /**< full-order constraint matrix (element dofs to local numbering) */
Mat Cp; /**< preconditioning constraint matrix (element dofs to local numbering, as sparse as possible) */
Vec weight; /**< Vector in global space, used to compensate for overcounting after local to global */
s_dFSWorkspace *workspace;
void *data;
};

PETSC_EXTERN_CXX_END
Expand Down
1 change: 1 addition & 0 deletions include/private/jacimpl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ struct _dRuleOps {
* nodes in each direction, weights in
* each direction. Does not copy, may not
* be implemented. */
dErr (*computeGeometry)(dRule,const dReal[restrict][3],dReal[restrict][3],dReal[restrict][3][3],dReal[restrict]);
};

/**
Expand Down
37 changes: 26 additions & 11 deletions src/fs/impls/cont/cont.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,10 @@ static dErr dFSBuildSpace_Cont(dFS fs)
/* MeshListEH v=MLZ,e=MLZ,f=MLZ,r=MLZ; */
/* MeshListInt in=MLZ,rfo=MLZ,feo=MLZ,rdegree=MLZ; */
/* dIInt nf; */
dInt i,j,gcnt,ghcnt,*inodes,*xnodes,*deg,*ghostind,nregions;
dInt *xind,*xstart,*istart,xcnt,*nnz,*pnnz;
dEntTopology *regTopo;
dInt i,j,gcnt,ghcnt,*inodes,*xnodes,*deg,*rdeg,*ghostind,nregions;
dInt *xind,*xstart,*istart,xcnt,*nnz,*pnnz,*regRDeg,*regBDeg;
dMeshEH *regEnts;
dEntStatus *status;
dErr err;

Expand All @@ -90,8 +92,9 @@ static dErr dFSBuildSpace_Cont(dFS fs)
err = dMeshGetAdjacency(mesh,fs->active,&ma);dCHK(err);
err = dFSContPropogateDegree(fs,&ma);dCHK(err);

err = dMallocA4(ma.nents*3,&deg,ma.nents,&inodes,ma.nents,&xnodes,ma.nents,&status);dCHK(err);
err = dMeshTagGetData(mesh,fs->degreetag,ma.ents,ma.nents,deg,3*ma.nents,dDATA_INT);
err = dMallocA5(ma.nents*3,&deg,ma.nents*3,&rdeg,ma.nents,&inodes,ma.nents,&xnodes,ma.nents,&status);dCHK(err);
err = dMeshTagGetData(mesh,fs->degreetag,ma.ents,ma.nents,deg,3*ma.nents,dDATA_INT);dCHK(err);
err = dMeshTagGetData(mesh,fs->ruletag,ma.ents,ma.nents,rdeg,3*ma.nents,dDATA_INT);dCHK(err);
/* Fill the arrays \a inodes and \a xnodes with the number of interior and expanded nodes for each
* (topology,degree) pair */
err = dJacobiGetNodeCount(fs->jacobi,ma.nents,ma.topo,deg,inodes,xnodes);dCHK(err);
Expand Down Expand Up @@ -169,11 +172,17 @@ static dErr dFSBuildSpace_Cont(dFS fs)
*/

nregions = ma.toff[dTYPE_REGION+1] - ma.toff[dTYPE_REGION];
err = dMallocA2(nregions,&xind,nregions+1,&xstart);dCHK(err);
err = dMallocA6(nregions,&xind,nregions+1,&xstart,nregions,&regTopo,nregions*3,&regRDeg,nregions*3,&regBDeg,nregions,&regEnts);dCHK(err);
for (i=0,xcnt=0; i<nregions; i++) {
const dInt ii = ma.toff[dTYPE_REGION] + i;
xind[i] = ii; /* Index of entity in full arrays */
xstart[i] = xcnt; /* first node on this entity */
regTopo[i] = ma.topo[ii];
regEnts[i] = ma.ents[ii];
for (j=0; j<3; j++) {
regRDeg[i*3+j] = dMaxInt(rdeg[ii*3+j],deg[ii*3+j]+2); /* Use a slightly stronger quadrature rule than required to be Hausdorff */
regBDeg[i*3+j] = deg[ii*3+j];
}
xcnt += xnodes[ii];
}
fs->m = xstart[nregions] = xcnt;
Expand All @@ -187,23 +196,29 @@ static dErr dFSBuildSpace_Cont(dFS fs)
err = dFree2(nnz,pnnz);dCHK(err);

err = dJacobiAddConstraints(fs->jacobi,nregions,xind,xstart,istart,deg,&ma,fs->C,fs->Cp);dCHK(err);
err = dFree2(xind,xstart);dCHK(err);
err = dFree(istart);dCHK(err);
err = dFree5(deg,rdeg,inodes,xnodes,status);dCHK(err);
err = dMeshRestoreAdjacency(mesh,fs->active,&ma);dCHK(err);

err = MatAssemblyBegin(fs->C,MAT_FINAL_ASSEMBLY);dCHK(err);
err = MatAssemblyBegin(fs->Cp,MAT_FINAL_ASSEMBLY);dCHK(err);
err = MatAssemblyEnd(fs->C,MAT_FINAL_ASSEMBLY);dCHK(err);
err = MatAssemblyEnd(fs->Cp,MAT_FINAL_ASSEMBLY);dCHK(err);

/* FIXME: get EFS and Rule */
/* Get Rule and EFS for domain ents. */
fs->nelem = nregions;
err = dMallocA3(nregions,&fs->rule,nregions,&fs->efs,nregions+1,&fs->off);dCHK(err); /* Will be freed by FS */
err = dMemcpy(fs->off,xstart,(nregions+1)*sizeof(xstart[0]));dCHK(err);
err = dJacobiGetRule(fs->jacobi,nregions,regTopo,regRDeg,fs->rule);dCHK(err);
err = dJacobiGetEFS(fs->jacobi,nregions,regTopo,regBDeg,fs->rule,fs->efs);dCHK(err);
err = dMeshGetVertexCoords(mesh,nregions,regEnts,&fs->vtxoff,&fs->vtx);dCHK(err); /* Should be restored by FS on destroy */
err = dFree6(xind,xstart,regTopo,regRDeg,regBDeg,regEnts);dCHK(err);

err = dFree(istart);dCHK(err);
err = dFree4(deg,inodes,xnodes,status);dCHK(err);
err = dMeshRestoreAdjacency(mesh,fs->active,&ma);dCHK(err);
dFunctionReturn(0);
}

/**
* Create the private structure used by a continuous galerkin function space.
* Create the private structure used by a continuous Galerkin function space.
*
* This function does not allocate the constraint matrices.
*
Expand Down
3 changes: 1 addition & 2 deletions src/fs/impls/cont/cont.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,7 @@ PETSC_EXTERN_CXX_BEGIN
extern dErr dFSCreate_Cont(dFS);

typedef struct {
dBool usecmatrix;
dMeshTag ruletag,degreetag;
dTruth usecmatrix;

/*
* Since we don't do matrix-free computations with these elements, they don't need to be persistant (i.e. defined
Expand Down
138 changes: 129 additions & 9 deletions src/fs/interface/fs.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,10 +66,9 @@ static inline dErr dBdyIntersect(dBdyType *a,dBdyType b)
dFunctionReturn(0);
}

/**
* Get an array of boundary types associated with a list of entities
*
* @param fs
/** Get an array of boundary types associated with a list of entities
*
* @param fs
* @param[in] nents Number of entities
* @param[in] ents handles
* @param[out] btype array of boundary types, should be at least \a nents in length
Expand Down Expand Up @@ -178,6 +177,13 @@ dErr dFSDestroy(dFS fs)
if (fs->ops->impldestroy) {
err = (*fs->ops->impldestroy)(fs);dCHK(err);
}
if (fs->workspace) {
s_dFSWorkspace *w = fs->workspace;
err = dFree7(w->q,w->jinv,w->jw,w->u,w->v,w->du,w->dv);dCHK(err);
err = dFree(fs->workspace);dCHK(err);
}
err = dFree3(fs->rule,fs->efs,fs->off);dCHK(err);
err = dMeshRestoreVertexCoords(fs->mesh,fs->nelem,NULL,&fs->vtxoff,&fs->vtx);dCHK(err);
err = PetscHeaderDestroy(fs);dCHK(err);
dFunctionReturn(0);
}
Expand All @@ -196,10 +202,7 @@ dErr dFSBuildSpace(dFS fs)

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
/*
* FIXME:
*
*/
if (fs->spacebuilt) dERROR(1,"The space is already built, rebuilding is not implemented");
if (fs->ops->buildspace) {
err = (*fs->ops->buildspace)(fs);dCHK(err);
}
Expand All @@ -214,7 +217,7 @@ dErr dFSBuildSpace(dFS fs)
err = VecGhostUpdateEnd(g,ADD_VALUES,SCATTER_FORWARD);dCHK(err);
err = VecDestroy(x);dCHK(err);

fs->spacebuilt = true;
fs->spacebuilt = dTRUE;
dFunctionReturn(0);
}

Expand Down Expand Up @@ -301,6 +304,123 @@ dErr dFSExpandedToGlobal(dFS fs,Vec x,InsertMode imode,Vec g)
dFunctionReturn(0);
}

/** Updates the global values, does \b not broadcast the global values back to the ghosts */
dErr dFSExpandedToGlobalBegin(dFS fs,Vec x,InsertMode imode,Vec g)
{
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidHeader(x,VEC_COOKIE,2);
dValidHeader(g,VEC_COOKIE,4);
err = dFSExpandedToGlobal(fs,x,imode,g);dCHK(err);
err = VecGhostUpdateBegin(g,imode,SCATTER_REVERSE);dCHK(err);
dFunctionReturn(0);
}

dErr dFSExpandedToGlobalEnd(dFS fs,Vec x,InsertMode imode,Vec g)
{
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
dValidHeader(x,VEC_COOKIE,2);
dValidHeader(g,VEC_COOKIE,4);
err = VecGhostUpdateEnd(g,imode,SCATTER_REVERSE);dCHK(err);
dFunctionReturn(0);
}

dErr dFSGetElements(dFS fs,dInt *n,dInt **off,s_dRule **rule,s_dEFS **efs,dInt **geomoff,dReal (**geom)[3])
{

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
if (n) *n = fs->nelem;
if (off) *off = fs->off;
if (rule) *rule = fs->rule;
if (efs) *efs = fs->efs;
if (geomoff) *geomoff = fs->vtxoff;
if (geom) *geom = fs->vtx;
dFunctionReturn(0);
}

dErr dFSRestoreElements(dFS fs,dInt *n,dInt **off,s_dRule **rule,s_dEFS **efs,dInt **geomoff,dReal (**geom)[3])
{

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
if (n) *n = 0;
if (off) *off = NULL;
if (rule) *rule = NULL;
if (efs) *efs = NULL;
if (geomoff) *geomoff = NULL;
if (geom) *geom = NULL;
dFunctionReturn(0);
}

/** Get arrays sufficiently large to hold the necessary quantities on the highest order element present in the mesh.
*
* @param fs
* @param q Pointer which will hold array of quadrature points, in physical space (not just the local tensor product)
* @param jinv Inverse of element Jacobian evaluated at quadrature points, normally just passed on to dEFSApply
* @param jw Array to store jacobian determinant times quadrature weights at quadrature points
* @param u first array to hold D values per quadrature point
* @param v second array to hold D values per quadrature point
* @param du first array to hold 3*D values per quadrature point
* @param dv second array to hold 3*D values per quadrature point
*/
dErr dFSGetWorkspace(dFS fs,dReal (**q)[3],dReal (**jinv)[3][3],dReal **jw,dScalar **u,dScalar **v,dScalar **du,dScalar **dv)
{
s_dFSWorkspace *w;
dErr err;

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
if (!fs->workspace) {
dInt Q = 0, D = fs->D;
err = dNew(s_dFSWorkspace,&w);dCHK(err);
for (dInt i=0; i<fs->nelem; i++) {
dInt nnodes;
err = dRuleGetSize(&fs->rule[i],NULL,&nnodes);dCHK(err);
if (nnodes > Q) Q = nnodes;
}
err = dMallocA7(Q,&w->q,Q,&w->jinv,Q,&w->jw,Q*D,&w->u,Q*D,&w->v,Q*D*3,&w->du,Q*D*3,&w->dv);dCHK(err);
fs->workspace = w;
}
w = fs->workspace;
if (q) *q = w->q;
if (jinv) *jinv = w->jinv;
if (jw) *jw = w->jw;
if (u) *u = w->u;
if (v) *v = w->v;
if (du) *du = w->du;
if (dv) *dv = w->dv;
dFunctionReturn(0);
}

/* These arrays are persistent for the life of the dFS so we just nullify the pointers */
dErr dFSRestoreWorkspace(dFS fs,dReal (**q)[3],dReal (**jinv)[3][3],dReal **jw,dScalar **u,dScalar **v,dScalar **du,dScalar **dv)
{

dFunctionBegin;
dValidHeader(fs,dFS_COOKIE,1);
#define dCHECK_NULLIFY(var) do { \
if (var) { \
if (*var == fs->workspace->var) *var = NULL; \
else dERROR(1,"Unable to restore array " #var " because it was not gotten or has been modified"); \
} \
} while (false)
dCHECK_NULLIFY(q);
dCHECK_NULLIFY(jinv);
dCHECK_NULLIFY(jw);
dCHECK_NULLIFY(u);
dCHECK_NULLIFY(v);
dCHECK_NULLIFY(du);
dCHECK_NULLIFY(dv);
#undef dCHECK_NULLIFY
dFunctionReturn(0);
}

dErr dFSGetMatrix(dFS fs,const MatType mtype,Mat *J)
{
dErr err;
Expand Down

0 comments on commit d022dba

Please sign in to comment.