Remove the constant pressure mode from computation of errors.

jedbrown · Apr 29, 2011 · 25fcbf7 · 25fcbf7
1 parent 73b6afc
commit 25fcbf7
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 3 deletions.
diff --git a/src/fs/tests/stokes.c b/src/fs/tests/stokes.c
@@ -819,6 +819,7 @@ static dErr StokesErrorNorms(Stokes stk,Vec gx,dReal errorNorms[3],dReal gerrorN
 {
   dErr             err;
   Vec              Coords,gxu,gxp;
+  PetscScalar      volume = 0,pressureshift = 0;
   dRulesetIterator iter;
 
   dFunctionBegin;
@@ -827,15 +828,34 @@ static dErr StokesErrorNorms(Stokes stk,Vec gx,dReal errorNorms[3],dReal gerrorN
   err = StokesExtractGlobalSplit(stk,gx,&gxu,&gxp);dCHK(err);
   err = StokesGetRegionIterator(stk,EVAL_FUNCTION,&iter);dCHK(err);
   err = dFSGetGeometryVectorExpanded(stk->fsu,&Coords);dCHK(err);
+  if (stk->alldirichlet) { // Do a volume integral of the exact solution to that we can remove the constant pressure mode
+    err = dRulesetIteratorStart(iter, Coords,dFS_INHOMOGENEOUS,NULL, NULL,NULL, NULL,NULL);dCHK(err);
+    while (dRulesetIteratorHasPatch(iter)) {
+      const dReal *jw;
+      const dScalar (*x)[3];
+      dInt Q;
+      err = dRulesetIteratorGetPatchApplied(iter,&Q,&jw, (dScalar**)&x,NULL,NULL,NULL, NULL,NULL,NULL,NULL, NULL,NULL,NULL,NULL);dCHK(err);
+      for (dInt i=0; i<Q; i++) {
+        dScalar uu[3],duu[9],pp[1],dpp[3];
+        stk->exact.solution(&stk->exactctx,x[i],uu,duu,pp,dpp);
+        volume += jw[i];
+        pressureshift += pp[0] * jw[i];
+      }
+      err = dRulesetIteratorNextPatch(iter);dCHK(err);
+    }
+    err = dRulesetIteratorFinish(iter);dCHK(err);
+    pressureshift /= volume;
+  }
   err = dRulesetIteratorStart(iter, Coords,dFS_INHOMOGENEOUS,NULL, gxu,dFS_INHOMOGENEOUS,NULL, gxp,dFS_INHOMOGENEOUS,NULL);dCHK(err);
   while (dRulesetIteratorHasPatch(iter)) {
     const dReal *jw;
     const dScalar (*x)[3],(*dx)[9],(*u)[3],(*du)[9],(*p)[1];
     dInt Q;
-    err = dRulesetIteratorGetPatchApplied(iter,&Q,&jw, (dScalar**)&x,(dScalar**)&dx,NULL,NULL, (const dScalar**)&u,(const dScalar**)&du,NULL,NULL, (const dScalar**)&p,NULL,NULL,NULL);dCHK(err);dCHK(err);
+    err = dRulesetIteratorGetPatchApplied(iter,&Q,&jw, (dScalar**)&x,(dScalar**)&dx,NULL,NULL, (const dScalar**)&u,(const dScalar**)&du,NULL,NULL, (const dScalar**)&p,NULL,NULL,NULL);dCHK(err);
     for (dInt i=0; i<Q; i++) {
       dScalar uu[3],duu[9],pp[1],dpp[3];
       stk->exact.solution(&stk->exactctx,x[i],uu,duu,pp,dpp);
+      pp[0] -= pressureshift;
       err = dNormsUpdate(errorNorms,gerrorNorms,jw[i],3,uu,u[i],duu,du[i]);dCHK(err);
       err = dNormsUpdate(perrorNorms,NULL,jw[i],1,pp,p[i],NULL,NULL);dCHK(err);
     }

diff --git a/src/fs/tests/stokesexact.py b/src/fs/tests/stokesexact.py
@@ -67,14 +67,14 @@ def solution(self, x,y,z, a,b,c):
                         1 * (c-1) * z,
                         0.25*(cos(pi*x) + cos(pi*y)) + 10*(x+y)])
 class StokesExact_1(StokesExact):
-    'From Larin & Reusken, 2009'
+    'From Larin & Reusken, 2009, with pressure shifted by a constant'
     def solution(self, x,y,z, a,b,c):
         from sympy import sin,cos,pi
         xx, yy, zz = (pi*s/2 for s in [x,y,z])
         return Matrix([+a/3 * sin(xx) * sin(yy) * sin(zz),
                        -b/3 * cos(xx) * cos(yy) * sin(zz),
                        -c*2/3 * cos(xx) * sin(yy) * cos(zz),
-                        cos(xx) * sin(yy) * sin(zz)])
+                        1 + cos(xx) * sin(yy) * sin(zz)])
 class StokesExact_2(StokesExact):
     def solution(self, x,y,z, a,b,c):
         return Matrix([a*z**3,