Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

Minor cleanup, fleshed out README.

  • Loading branch information...
commit fcfe7a61da774f900cc5cefa2820944f9d0d812c 1 parent acd5393
@jedbrown authored
Showing with 72 additions and 20 deletions.
  1. +2 −1  Makefile
  2. +54 −1 README
  3. +1 −1  elliptic.C
  4. +15 −17 stokes.C
View
3  Makefile
@@ -22,7 +22,8 @@ poisson : poisson.o chebyshev.o chkopts
${CLINKER} -o poisson poisson.o chebyshev.o -lfftw3 ${PETSC_LIB}
elliptic : elliptic.o chebyshev.o chkopts
- ${CLINKER} -o elliptic elliptic.o chebyshev.o ${PROF} -lfftw3 ${PETSC_LIB}
+ @-echo 'Linking elliptic'
+ @-${CLINKER} -o elliptic elliptic.o chebyshev.o ${PROF} -lfftw3 ${PETSC_LIB}
stokes : stokes.o chebyshev.o chkopts
@-echo 'Linking stokes'
View
55 README
@@ -1,7 +1,9 @@
Experiments with matrix-free Chebyshev collocation for nonlinear elliptic/Stokes
problems with fully iterative solutions.
-You will need FFTW and a recent PETSc to build.
+You will need FFTW and a recent PETSc to build. If you are not using petsc-dev,
+switch the include line in the Makefile. If you do not have CppAD, define
+WITH_CPPAD to be 0 in stokes.C.
The numerical method is spectral collocation on the Chebyshev-Gauss-Lobatto
nodes. Derivatives are applied by DCT using FFTW. For preconditioning, a low
@@ -13,3 +15,54 @@ number of required iterations is nearly independent of polynomial order. Thus,
for a 3D problem with polynomial order P in each cartesion direction, there are
P^3 nodes and O(P^3 log P) operations required when using multigrid
preconditioning. This is optimal for a spectral method.
+
+The scalar elliptic implementation is truly arbitrary dimensional. For instance
+
+ ./elliptic -dim 12,12,12,12,12 -pc_type hypre -exact 2 -ksp_monitor -ksp_rtol 1e-10
+
+solves the Poisson problem in 5 dimensions with multigrid preconditioning,
+comparing to an exact solution with inhomogenous Dirichlet data.
+
+For the Stokes problem, the we must precondition the saddle point problem. This
+is done by approximately solving with the Schur complement (a block LU
+decomposition). There are KSP objects to control. The top level will typically
+need to be FGMRES since the preconditioner is usually not linear. The Schur
+complement is controlled by -schur_ksp_xxx. For slow flow problems, the Schur
+complement is usually preconditioned with the mass matrix. For a collocation
+method, the mass matrix is the identity, so we do not precondition it. Perhaps
+we can do something clever for the power law rheology which increases the
+condition number. Part of the block LU involves a solve with the (viscous)
+velocity problem. This is controlled by -vel_xx. The preconditioner matrix is
+the sparse approximation. Typically, we will use AMG and a few iterations (~5)
+of GMRES on this problem. Applying the Schur complement also involves an
+(approximate) solve with the velocity system. This solve is controlled by
+-svel_xx. Here, it is less beneficial to do a good solve, so we may do just one
+V-cycle on the sparse approximation. Putting this all together, we have
+something like:
+
+ ./stokes -exact 2 -cont0 1 -schur_ksp_max_it 3 -vel_ksp_max_it 4 -vel_pc_type hypre -svel_ksp_type preonly -svel_pc_type hypre -ksp_type fgmres -ksp_monitor -dim 20,20,20 -ksp_rtol 1e-10
+
+When the power law rheology has very weak regularization, we must do a
+continuation for the Newton iteration to converge. The option -cont0 <n> gives
+the starting value of the continuation (usually 0 to solve the linear problem)
+and -cont <n> gives the final value. Use -rheology 1 to use power law rheology.
+Note that the error for exact solutions no longer applies since they are for
+linear viscosity.
+
+ ./stokes -exact 2 -cont 4 -rheology 1 -eps 1e-4 -exponent 3 -schur_ksp_max_it 3 -vel_ksp_max_it 4 -vel_pc_type hypre -svel_ksp_type preonly -svel_pc_type hypre -ksp_type fgmres -ksp_monitor -snes_monitor -dim 20,20,20
+
+
+For the preconditioning matrix, there are a few options, although the simple
+finite difference poisson-type preconditioner seems pretty good. I suspect the
+Q1 finite element version will do better with a high order Galerkin method.
+Interestingly, the finite difference via automatic differentiation does not
+typically work better than the simple version. Also, there seems to be a bug in
+CppAD which makes it break on certain problems. Note that just subsampling the
+spectral matrix with specified stencil provides a very poor preconditioner.
+Choose velocity preconditioning matrices with -pcvel <n>.
+
+The current approach for Neumann and mixed boundary conditions is
+broken/incomplete. While the Neumann condition works, it destroys the
+conditioning of the system so convergence is slow. I am not sure that the mixed
+condition is even correct and convergence is terrible. A suitable outflow
+boundary is also needed, but this is not implemented.
View
2  elliptic.C
@@ -140,7 +140,7 @@ int main(int argc,char **args)
ac->dim[0] = 8; ac->dim[1] = 6; // default dimension extent
ac->debug = 0; ac->exact = 0; ac->gamma = 0.0; ac->exponent = 2.0;
ierr = PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elliptic problem options", ""); CHKERRQ(ierr);
- ierr = PetscOptionsIntArray("-dims", "list of dimension extent", "elliptic.C", ac->dim, &ac->d, &flag); CHKERRQ(ierr);
+ ierr = PetscOptionsIntArray("-dim", "list of dimension extent", "elliptic.C", ac->dim, &ac->d, &flag); CHKERRQ(ierr);
if (!flag) ac->d = 2;
ierr = PetscOptionsInt("-debug", "debugging level", "elliptic.C", ac->debug, &ac->debug, PETSC_NULL); CHKERRQ(ierr);
ierr = PetscOptionsInt("-exact", "exact solution type", "elliptic.C", ac->exact, &ac->exact, PETSC_NULL); CHKERRQ(ierr);
View
32 stokes.C
@@ -1,7 +1,5 @@
static char help[] = "Stokes problem with non-Newtonian rheology via Chebyshev collocation.\n";
-//#define InFunction printf("%s\n", __FUNCT__)
-#define InFunction
#define SOLVE 1
#define CHECK_EXACT 1
#define WITH_CPPAD 1
@@ -161,11 +159,11 @@ int main(int argc,char **args)
ierr = PetscOptionsGetInt(PETSC_NULL, "-pcvel", &t, &flag);CHKERRQ(ierr);
if (!flag) t = 0;
switch (t) {
- case 0: ierr = PCShellSetSetUp(pc, StokesPCSetUp0);CHKERRQ(ierr); break;
- case 1: ierr = PCShellSetSetUp(pc, StokesPCSetUp1);CHKERRQ(ierr); break;
- case 2: ierr = PCShellSetSetUp(pc, StokesPCSetUp2);CHKERRQ(ierr); break;
+ case 0: ierr = PCShellSetSetUp(pc, StokesPCSetUp0);CHKERRQ(ierr); break; // Simple finite difference
+ case 1: ierr = PCShellSetSetUp(pc, StokesPCSetUp1);CHKERRQ(ierr); break; // Finite element
+ case 2: ierr = PCShellSetSetUp(pc, StokesPCSetUp2);CHKERRQ(ierr); break; // Subsampling the spectral matrix using coloring
#if (WITH_CPPAD)
- case 3: ierr = PCShellSetSetUp(pc, StokesPCSetUp3);CHKERRQ(ierr); break;
+ case 3: ierr = PCShellSetSetUp(pc, StokesPCSetUp3);CHKERRQ(ierr); break; // Finite difference with automatic differentiation
#endif
default: SETERRQ1(1, "pcvel type number %d not implemented", t);CHKERRQ(ierr);
}
@@ -640,7 +638,7 @@ PetscErrorCode StokesFunction(SNES snes, Vec xG, Vec yG, void *ctx)
PetscReal **v, **strain, *eta, *deta;
PetscErrorCode ierr;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
d = c->numDims; n = productInt(d, c->dim);
xL = c->workV[0]; yL = c->workV[1]; V = &c->workV[2]; W = &c->workV[2+d];
ierr = VecScatterBegin(c->scatterGP, xG, c->pG0, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr);
@@ -1782,7 +1780,7 @@ PetscErrorCode StokesVecView(Vec v, PetscInt nodes, PetscInt pernode, PetscInt p
PetscErrorCode StokesRheologyLinear(PetscInt d, PetscReal gamma, PetscReal *eta, PetscReal *deta, void *ctx)
{
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
*eta = 1.0; *deta = 0.0;
PetscFunctionReturn(0);
}
@@ -1810,7 +1808,7 @@ PetscErrorCode StokesRheologyPower(PetscInt d, PetscReal gamma, PetscReal *eta,
PetscErrorCode StokesExact0(PetscInt d, PetscReal *coord, PetscReal *value, PetscReal *rhs, void *ctx)
{
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
if (value) {
for (int i=0; i < d+1; i++) value[i] = 0.0;
}
@@ -1827,7 +1825,7 @@ PetscErrorCode StokesExact1(PetscInt d, PetscReal *coord, PetscReal *value, Pets
const PetscReal eta = 1.0;
PetscReal u, v, p;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
if (d > 3) SETERRQ2(1, "%s only implemented for dimension 2 and 3 but %d given", __FUNCT__, d);
u = sin(0.5 * PETSC_PI * coord[0]) * cos(0.5 * PETSC_PI * coord[1]);
v = -cos(0.5 * PETSC_PI * coord[0]) * sin(0.5 * PETSC_PI * coord[1]);
@@ -1855,7 +1853,7 @@ PetscErrorCode StokesExact2(PetscInt d, PetscReal *coord, PetscReal *value, Pets
const PetscReal eta = 1.0;
PetscReal u, v, p;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
if (d > 3) SETERRQ2(1, "%s only implemented for dimension 2 but %d given", __FUNCT__, d);
u = sin(0.5 * PETSC_PI * coord[0]) * cos(0.5 * PETSC_PI * coord[1]);
v = -cos(0.5 * PETSC_PI * coord[0]) * sin(0.5 * PETSC_PI * coord[1]);
@@ -1880,7 +1878,7 @@ PetscErrorCode StokesExact3(PetscInt d, PetscReal *coord, PetscReal *value, Pets
const PetscReal eta = 1.0;
PetscReal u, v, p;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
if (d != 2) SETERRQ2(1, "%s only implemented for dimension 2 but %d given", __FUNCT__, d);
u = coord[1] + 1.0;
v = 0.0;
@@ -1904,7 +1902,7 @@ PetscErrorCode StokesDirichlet(PetscInt d, PetscReal *coord, PetscReal *normal,
StokesExactBoundaryCtx *ctx = (StokesExactBoundaryCtx *)void_ctx;
PetscErrorCode ierr;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
*type = DIRICHLET;
ierr = ctx->exact(d, coord, value, PETSC_NULL, ctx->exactCtx);CHKERRQ(ierr);
PetscFunctionReturn(0);
@@ -1919,7 +1917,7 @@ PetscErrorCode StokesBoundary1(PetscInt d, PetscReal *coord, PetscReal *normal,
bool inside;
PetscErrorCode ierr;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
inside = false;
for (int i=0; i < d-1; i++) inside = inside || (PetscAbs(coord[i]) < 0.999);
if (coord[d-1] > 0.999 && inside) { // Impose condition at the 'surface'
@@ -1960,7 +1958,7 @@ PetscErrorCode StokesBoundary2(PetscInt d, PetscReal *coord, PetscReal *normal,
bool inside;
PetscErrorCode ierr;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
inside = false;
for (int i=0; i < d-1; i++) inside = inside || (PetscAbs(coord[i]) < 0.999);
if (coord[d-1] > 0.999 && inside) { // Impose condition at the 'surface'
@@ -2003,7 +2001,7 @@ PetscErrorCode StokesBoundary3(PetscInt d, PetscReal *coord, PetscReal *normal,
bool inside;
PetscErrorCode ierr;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
inside = false;
for (int i=0; i < d-1; i++) inside = inside || (PetscAbs(coord[i]) < 0.999);
inside = true;
@@ -2033,7 +2031,7 @@ PetscErrorCode StokesBoundary4(PetscInt d, PetscReal *coord, PetscReal *normal,
bool inside;
PetscErrorCode ierr;
- PetscFunctionBegin; InFunction;
+ PetscFunctionBegin;
*type = DIRICHLET;
for (int i=0; i < d; i++) value[i] = 0.0;
if (coord[d-2] < -0.999) value[d-2] = 1 - 0.25 * PetscSqr(coord[d-1]-1);
Please sign in to comment.
Something went wrong with that request. Please try again.