Fix jacobi tests for dEFSSparseSetUp_Tensor

jedbrown · Dec 22, 2010 · 94b220f · 94b220f
1 parent d150390
commit 94b220f
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 13 deletions.
diff --git a/src/fs/tests/CMakeLists.txt b/src/fs/tests/CMakeLists.txt
@@ -7,7 +7,7 @@ dohp_link_executable (cunit cunit.c)
 dohp_link_executable (elast elast.c)
 dohp_link_executable (stokes stokes.c)
 
-dohp_add_test (fs-mf-projection 1 fs-ex1 -snes_mf -ksp_type minres -const_BDeg 6 -snes_monitor -ksp_converged_reason -dfs_ordering_type natural)
+dohp_add_test (fs-mf-projection 1 fs-ex1 -snes_mf -ksp_type minres -const_BDeg 5 -snes_monitor_short -ksp_converged_reason -dfs_ordering_type natural)
 dohp_add_test (fs-mf-op-proj-0  1 fs-ex1 -snes_mf_operator -ksp_type minres -pc_type jacobi -const_BDeg 10 -snes_monitor -ksp_converged_reason -ksp_rtol 1e-10 -frequency 5,4,3 -require_ptwise 2e-6,3e-6,2.3e-5 -require_grad 7e-7,3e-6,2e-4 -dfs_ordering_type natural)
 dohp_add_test (fs-mf-op-proj-1  1 fs-ex1 -const_BDeg 7 -snes_mf_operator -ksp_converged_reason -pc_type jacobi -ksp_type minres -ksp_rtol 1e-10 -snes_monitor -frequency 2,2,2 -require_ptwise 2e-6,2e-6,9e-6 -require_grad 3e-6,7e-6,1.1e-4 -dfs_ordering_type natural)
 

diff --git a/src/fs/tests/refout/fs-mf-projection.refout b/src/fs/tests/refout/fs-mf-projection.refout
@@ -1,12 +1,9 @@
 expanded size 432, global size 396, range 0 .. 396 
 |x| = 432.000000, |g| = 432.000000, |y| = 504.000000
-  0 SNES Function norm 3.536327126455e-01 
-Linear solve converged due to CONVERGED_RTOL iterations 96
-  1 SNES Function norm 3.505422504790e-06 
-Linear solve converged due to CONVERGED_RTOL iterations 206
-  2 SNES Function norm 3.465700777016e-11 
-Linear solve converged due to CONVERGED_RTOL iterations 209
-  3 SNES Function norm 3.669444474224e-16 
-Algebraic projection residual  |x|_1 4.37e-15  |x|_2 3.67e-16  |x|_inf 1.27e-16
-Pointwise projection residual  |x|_1 2.79e-07  |x|_2 3.77e-07  |x|_inf 1.88e-06
-Gradient ptwise proj residual  |x|_1 6.07e-08  |x|_2 2.00e-07  |x|_inf 3.28e-06
+  0 SNES Function norm 0.353633 
+  Linear solve converged due to CONVERGED_RTOL iterations 97
+  1 SNES Function norm 3.22358e-06 
+  Linear solve converged due to CONVERGED_RTOL iterations 206
+  2 SNES Function norm 3.263e-11 
+  Linear solve converged due to CONVERGED_RTOL iterations 209
+  3 SNES Function norm < 1.e-11
diff --git a/src/jacobi/impls/tensor/tensor.c b/src/jacobi/impls/tensor/tensor.c
@@ -495,6 +495,24 @@ static dErr TensorGetBasis(dJacobi_Tensor *tnsr,dInt rsize,const dReal rcoord[],
   dFunctionReturn(0);
 }
 
+static dErr dEFSSparseSetUp_Tensor_Default(dEFS_Tensor *tefs)
+{
+  struct dEFS_TensorSparse *sparse = &tefs->sparse;
+  const dReal *interp,*deriv;
+  dInt Q,P;
+  dErr err;
+
+  dFunctionBegin;
+  err = dEFSGetExplicit((dEFS)tefs,NULL,&Q,&P,&interp,&deriv);dCHK(err);
+  sparse->npieces = 1;
+  err = dMallocA6(1,&sparse->Q,1,&sparse->P,1,&sparse->qidx,1,&sparse->eidx,1,&sparse->interp,1,&sparse->deriv);dCHK(err);
+  err = dMallocA4(Q,&sparse->qidx[0],P,&sparse->eidx[0],Q*P,&sparse->interp[0],Q*P*3,&sparse->deriv[0]);dCHK(err);
+  err = dMemcpy(sparse->interp[0],interp,Q*P*sizeof(interp[0]));dCHK(err);
+  err = dMemcpy(sparse->deriv[0],deriv,Q*P*3*sizeof(deriv[0]));dCHK(err);
+  err = dEFSRestoreExplicit((dEFS)tefs,NULL,&Q,&P,&interp,&deriv);dCHK(err);
+  dFunctionReturn(0);
+}
+
 static dErr dEFSSparseSetUp_Tensor(dEFS_Tensor *tefs)
 {
   struct dEFS_TensorSparse *sparse = &tefs->sparse;
@@ -506,8 +524,11 @@ static dErr dEFSSparseSetUp_Tensor(dEFS_Tensor *tefs)
     const dInt P[3] = {tefs->basis[0]->P,tefs->basis[1]->P,tefs->basis[2]->P},Q[3] = {tefs->basis[0]->Q,tefs->basis[1]->Q,tefs->basis[2]->Q};
     dInt N,QQ,qperpiece,eperpiece,*qidx,*eidx;
     dReal *interp,*deriv;
-    if (Q[0] != 2*P[0] || Q[1] != 2*P[1] || Q[2] != 2*P[2])
-      dERROR(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Sparse assembly not possible with this order-quadrature");
+    if (Q[0] != 2*P[0] || Q[1] != 2*P[1] || Q[2] != 2*P[2]) {
+      /* Sparse assembly not possible with this order quadrature */
+      err = dEFSSparseSetUp_Tensor_Default(tefs);dCHK(err);
+      dFunctionReturn(0);
+    }
     sparse->npieces = N = P[0] * P[1] * P[2];
     QQ = Q[0]*Q[1]*Q[2];
     qperpiece = 8;              /* Sloppy, assumes 2^3 quadrature on each piece */
@@ -562,6 +583,7 @@ static dErr dJacobiGetEFS_Tensor(dJacobi jac,dInt n,const dEntTopology topo[],co
           newefs->ops = *tnsr->efsOpsQuad;
           err = TensorGetBasis(tnsr,rsize[0],rcoord[0],dPolynomialOrder1D(order[i],0),&newefs->basis[0]);dCHK(err);
           err = TensorGetBasis(tnsr,rsize[1],rcoord[1],dPolynomialOrder1D(order[i],1),&newefs->basis[1]);dCHK(err);
+          err = dEFSSparseSetUp_Tensor(newefs);dCHK(err);
           break;
         case dTOPO_HEX:
           if (rdim != 3) dERROR(PETSC_COMM_SELF,1,"Incompatible Rule dim %d, expected 3",rdim);