Skip to content

Commit

Permalink
Stokes: add -gravity, project Dirichlet part
Browse files Browse the repository at this point in the history
  • Loading branch information
@jedbrown
jedbrown committed Apr 30, 2011
1 parent 0ccd02d commit f9fa034
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 21 deletions.
64 changes: 44 additions & 20 deletions src/fs/tests/stokes.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,12 +22,27 @@ static const char help[] = "Solve non-Newtonian Stokes problem using dual order

#include <stokesexact.h> /* Generated by stokesexact.py */

static void StokesExact_Gravity_Solution(const struct StokesExactCtx dUNUSED *exc,const dReal dUNUSED x[3],dScalar u[],dScalar du[],dScalar *p,dScalar dp[])
{ /* Defines inhomogeneous Dirichlet boundary conditions */
u[0] = u[1] = u[2] = 0;
for (dInt i=0; i<9; i++) du[i] = 0;
*p = 0;
for (dInt i=0; i<3; i++) dp[i] = 0;
}
static void StokesExact_Gravity_Forcing(const struct StokesExactCtx dUNUSED *exc,const struct StokesRheology *rheo,const dReal dUNUSED x[3],dScalar fu[],dScalar *fp)
{
fu[0] = 0;
fu[1] = 0;
fu[2] = -rheo->gravity;
fp[0] = 0;
}

#define ALEN(a) ((dInt)(sizeof(a)/sizeof(a)[0]))
static PetscLogEvent LOG_StokesShellMult;
typedef struct _p_Stokes *Stokes;
typedef struct _n_Stokes *Stokes;

struct StokesExact {
void (*solution)(const struct StokesExactCtx*,const dReal x[3],dScalar u[],dScalar *p,dScalar du[],dScalar dp[]);
void (*solution)(const struct StokesExactCtx*,const dReal x[3],dScalar u[],dScalar du[],dScalar *p,dScalar dp[]);
void (*forcing)(const struct StokesExactCtx*,const struct StokesRheology*,const dReal x[3],dScalar fu[],dScalar *fp);
};

Expand All @@ -50,19 +65,20 @@ static dErr StokesShellMatMultAdd_Bt(Mat A,Vec gx,Vec gy,Vec gz) {return StokesS
static dErr StokesShellMatMultAdd_B(Mat A,Vec gx,Vec gy,Vec gz) {return StokesShellMatMult_All_IorA(A,gx,gy,gz,ADD_VALUES,STOKES_MULT_B);}
static dErr MatGetVecs_Stokes(Mat,Vec*,Vec*);

struct _p_Stokes {
struct _n_Stokes {
MPI_Comm comm;
struct StokesRheology rheo;
struct StokesExact exact;
struct StokesExactCtx exactctx;
dBool reality;
dFS fsu,fsp;
Vec xu,xp,yu,yp;
Vec gvelocity,gpressure,gpacked;
IS ublock,pblock; /* Global index sets for each block */
IS lublock,lpblock; /* Local index sets for each block */
VecScatter extractVelocity,extractPressure;
dInt constBDeg,pressureCodim;
dBool errorview,cardinalMass;
dBool cardinalMass;
char mattype_A[256],mattype_D[256];
dQuadratureMethod function_qmethod,jacobian_qmethod;
dRulesetIterator regioniter[EVAL_UB];
Expand All @@ -77,7 +93,7 @@ static dErr StokesCreate(MPI_Comm comm,Stokes *stokes)

dFunctionBegin;
*stokes = 0;
err = dNew(struct _p_Stokes,&stk);dCHK(err);
err = dNew(struct _n_Stokes,&stk);dCHK(err);
stk->comm = comm;

stk->constBDeg = 3;
Expand Down Expand Up @@ -145,17 +161,19 @@ static dErr StokesSetFromOptions(Stokes stk)
err = PetscOptionsInt("-const_bdeg","Use constant isotropic degree on all elements","",stk->constBDeg,&stk->constBDeg,NULL);dCHK(err);
err = PetscOptionsInt("-pressure_codim","Reduce pressure space by this factor","",stk->pressureCodim,&stk->pressureCodim,NULL);dCHK(err);
err = PetscOptionsBool("-cardinal_mass","Assemble diagonal mass matrix","",stk->cardinalMass,&stk->cardinalMass,NULL);dCHK(err);
err = PetscOptionsBool("-error_view","View errors","",stk->errorview,&stk->errorview,NULL);dCHK(err);
err = PetscOptionsReal("-rheo_A","Rate factor (rheology)","",rheo->A,&rheo->A,NULL);dCHK(err);
err = PetscOptionsReal("-rheo_eps","Regularization (rheology)","",rheo->eps,&rheo->eps,NULL);dCHK(err);
err = PetscOptionsReal("-rheo_p","Power p=1+1/n where n is Glen exponent","",rheo->p,&rheo->p,NULL);dCHK(err);
err = PetscOptionsInt("-exact","Exact solution choice","",exact,&exact,NULL);dCHK(err);
err = PetscOptionsReal("-exact_a","First scale parameter","",exc->a,&exc->a,NULL);dCHK(err);
err = PetscOptionsReal("-exact_b","Second scale parameter","",exc->b,&exc->b,NULL);dCHK(err);
err = PetscOptionsReal("-exact_c","Third scale parameter","",exc->c,&exc->c,NULL);dCHK(err);
err = PetscOptionsReal("-exact_scale","Overall scale parameter","",exc->scale,&exc->scale,NULL);dCHK(err);
err = PetscOptionsList("-stokes_A_mat_type","Matrix type for velocity operator","",MatList,stk->mattype_A,stk->mattype_A,sizeof(stk->mattype_A),NULL);dCHK(err);
err = PetscOptionsList("-stokes_D_mat_type","Matrix type for velocity operator","",MatList,stk->mattype_D,stk->mattype_D,sizeof(stk->mattype_D),NULL);dCHK(err);
err = PetscOptionsReal("-gravity","Force with gravity equal to this value (non-dimensional density=1)","",rheo->gravity,&rheo->gravity,&stk->reality);dCHK(err);
if (!stk->reality) {
err = PetscOptionsInt("-exact","Exact solution choice","",exact,&exact,NULL);dCHK(err);
err = PetscOptionsReal("-exact_a","First scale parameter","",exc->a,&exc->a,NULL);dCHK(err);
err = PetscOptionsReal("-exact_b","Second scale parameter","",exc->b,&exc->b,NULL);dCHK(err);
err = PetscOptionsReal("-exact_c","Third scale parameter","",exc->c,&exc->c,NULL);dCHK(err);
err = PetscOptionsReal("-exact_scale","Overall scale parameter","",exc->scale,&exc->scale,NULL);dCHK(err);
}
err = PetscOptionsList("-Ap_mat_type","Matrix type for velocity operator","",MatList,stk->mattype_A,stk->mattype_A,sizeof(stk->mattype_A),NULL);dCHK(err);
err = PetscOptionsList("-Dp_mat_type","Matrix type for velocity operator","",MatList,stk->mattype_D,stk->mattype_D,sizeof(stk->mattype_D),NULL);dCHK(err);
err = PetscOptionsEnum("-stokes_f_qmethod","Quadrature method for residual evaluation/matrix-free","",dQuadratureMethods,(PetscEnum)stk->function_qmethod,(PetscEnum*)&stk->function_qmethod,NULL);dCHK(err);
err = PetscOptionsEnum("-stokes_jac_qmethod","Quadrature to use for Jacobian assembly","",dQuadratureMethods,(PetscEnum)stk->jacobian_qmethod,(PetscEnum*)&stk->jacobian_qmethod,NULL);dCHK(err);
{
Expand All @@ -168,7 +186,10 @@ static dErr StokesSetFromOptions(Stokes stk)
}
} err = PetscOptionsEnd();dCHK(err);

switch (exact) {
if (stk->reality) {
stk->exact.solution = StokesExact_Gravity_Solution;
stk->exact.forcing = StokesExact_Gravity_Forcing;
} else switch (exact) {
case 0:
stk->exact.solution = StokesExact_0_Solution;
stk->exact.forcing = StokesExact_0_Forcing;
Expand Down Expand Up @@ -1044,9 +1065,9 @@ int main(int argc,char *argv[])
MPI_Comm comm;
Mat J,Jp;
MatFDColoring fdcolor = NULL;
Vec r,x,soln;
Vec r,x,soln = NULL;
SNES snes;
dBool nocheck,check_null,compute_explicit,use_jblock,viewdhm;
dBool check_error,check_null,compute_explicit,use_jblock,viewdhm;
dErr err;

err = dInitialize(&argc,&argv,NULL,help);dCHK(err);
Expand All @@ -1060,7 +1081,8 @@ int main(int argc,char *argv[])
err = VecDuplicate(r,&x);dCHK(err);

err = PetscOptionsBegin(stk->comm,NULL,"Stokes solver options",__FILE__);dCHK(err); {
err = PetscOptionsBool("-nocheck_error","Do not compute errors","",nocheck=dFALSE,&nocheck,NULL);dCHK(err);
check_error = stk->reality ? dFALSE : dTRUE;
err = PetscOptionsBool("-check_error","Compute errors","",check_error,&check_error,NULL);dCHK(err);
err = PetscOptionsBool("-use_jblock","Use blocks to apply Jacobian instead of unified (more efficient) version","",use_jblock=dFALSE,&use_jblock,NULL);dCHK(err);
err = PetscOptionsBool("-viewdhm","View the solution","",viewdhm=dFALSE,&viewdhm,NULL);dCHK(err);
err = PetscOptionsBool("-check_null","Check that constant pressure really is in the null space","",check_null=dFALSE,&check_null,NULL);dCHK(err);
Expand Down Expand Up @@ -1099,7 +1121,7 @@ int main(int argc,char *argv[])
err = PCFieldSplitSetIS(pc,"p",stk->pblock);dCHK(err);
}
err = StokesGetSolutionVector(stk,&soln);dCHK(err);
{
if (!stk->reality) {
dReal nrm;
MatStructure mstruct;
Vec b;
Expand All @@ -1123,7 +1145,7 @@ int main(int argc,char *argv[])
err = StokesGetNullSpace(stk,&matnull);dCHK(err);
err = SNESGetKSP(snes,&ksp);dCHK(err);
err = KSPSetNullSpace(ksp,matnull);dCHK(err);
err = MatNullSpaceRemove(matnull,soln,NULL);dCHK(err);
if (soln) {err = MatNullSpaceRemove(matnull,soln,NULL);dCHK(err);}
err = MatNullSpaceDestroy(&matnull);dCHK(err);
}
if (check_null) {
Expand All @@ -1139,7 +1161,7 @@ int main(int argc,char *argv[])
err = KSPGetNullSpace(ksp,&matnull);dCHK(err); /* does not reference */
err = MatNullSpaceRemove(matnull,x,NULL);dCHK(err);
}
if (!nocheck) {
if (check_error) {
dReal anorm[3],inorm[3],enorm[3],gnorm[3],epnorm[3];
err = StokesErrorNorms(stk,x,enorm,gnorm,epnorm);dCHK(err);
err = dNormsAlgebraicScaled(anorm,r);dCHK(err);
Expand All @@ -1159,6 +1181,8 @@ int main(int argc,char *argv[])
err = PetscViewerFileSetName(view,"stokes.dhm");dCHK(err);
err = PetscViewerFileSetMode(view,FILE_MODE_WRITE);dCHK(err);
err = StokesExtractGlobalSplit(stk,x,&Xu,&Xp);dCHK(err);
err = dFSDirichletProject(stk->fsu,Xu,dFS_INHOMOGENEOUS);dCHK(err);
err = dFSDirichletProject(stk->fsp,Xp,dFS_INHOMOGENEOUS);dCHK(err);
err = VecView(Xu,view);dCHK(err);
err = VecView(Xp,view);dCHK(err);
err = PetscViewerDestroy(&view);dCHK(err);
Expand Down
2 changes: 1 addition & 1 deletion src/fs/tests/stokesexact.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ def solution(self, x,y,z, a,b,c):
solutions = [StokesExact_0(), StokesExact_1(), StokesExact_2(), StokesExact_3()]
with open('stokesexact.h', 'w') as fheader, open('stokesexact.c', 'w') as fimpl:
fheader.write('#include <dohptype.h>\n\n')
fheader.write('struct StokesRheology {dReal A,eps,p;};\n')
fheader.write('struct StokesRheology {dReal A,eps,p,gravity;};\n') # Gravity is not actually used by exact solutions
fheader.write('struct StokesExactCtx {dReal a,b,c,scale;};\n')
fimpl.write('\n'.join(['#include "stokesexact.h"', '#include <math.h>', '#ifndef M_PI', '# define M_PI 3.14159265358979323846', '#endif', '\n']))
for sol in solutions:
Expand Down

0 comments on commit f9fa034

Please sign in to comment.