Permalink
Browse files

Add plain C support, testing conditioning of perturbed matrices

Makefile: Allow a plain C version of PETSc to be used
stokes.C: Testing conditioning of perturbed matrices

Add option to filter matrix for testing effect of imposing
Dirichlet conditions by zeroing diagonal.  The option are only effective if
-pc_type 1 (Q1 finite element) preconditioning is chosen.  There are two new
options:

  -zeroN Int      number of rows and columns to eliminate
  -zeroV Scalar   value to place on the diagonal

We can test the conditioning of the velocity matrix with

  -vel_ksp_type preonly -vel_pc_type ksp -vel_ksp_ksp_rtol 1e-10 \
  -vel_ksp_ksp_gmres_restart 500 -vel_ksp_pc_type none \
  -vel_ksp_ksp_monitor_singular_value

and the conditioning of the Schur complement with

  -schur_pc_type none -svel_ksp_type preonly -svel_pc_type ksp \
  -svel_ksp_rtol 1e-10 -svel_ksp_pc_type ilu -schur_ksp_gmres_restart 500 \
  -schur_ksp_monitor_singular_value

Note that the divergence matrices are still applied to spectral order and
matrix-free.  Hence the procedure here is only a first-order approximation
to the effect on conditioning since we are not removing the appropriate rows
and columns from these matrix-free operators.
  • Loading branch information...
1 parent 7b5a751 commit c4ec55acd7bd70d2e532247daba08041d99be99a @jedbrown committed Aug 15, 2008
Showing with 20 additions and 3 deletions.
  1. +1 −0 Makefile
  2. +19 −3 stokes.C
View
@@ -5,6 +5,7 @@ CFLAGS+= ${PROF}
ETAGS= etags.emacs
include ${PETSC_DIR}/conf/base
#include ${PETSC_DIR}/bmake/common/base
+CLINKER= mpicxx
nk : nk.o chkopts
${CLINKER} -o nk nk.o ${PETSC_LIB}
View
@@ -29,8 +29,8 @@ typedef struct { // For higher dimensions, use a different value of `3'
} StokesBoundaryMixed;
typedef struct {
- PetscInt debug, cont0, cont;
- PetscReal hardness, exponent, regularization, gamma0, scaleM, scaleN, reynolds;
+ PetscInt debug, cont0, cont, zeroN;
+ PetscReal hardness, exponent, regularization, gamma0, scaleM, scaleN, reynolds, zeroV;
Rheology rheology;
ExactSolution exact;
BdyFunc boundary;
@@ -393,7 +393,7 @@ PetscErrorCode StokesProcessOptions(StokesCtx *ctx)
ierr = PetscMalloc(ctx->numDims*sizeof(PetscInt), &ctx->dim);CHKERRQ(ierr);
exact = 0; boundary = 0; rheology = 0;
opt->debug = 0; opt->hardness = 1.0; opt->exponent = 1.0; opt->regularization = 1.0; opt->gamma0 = 1.0; opt->cont0 = 0; opt->cont = 1;
- opt->scaleM = 1.0; opt->scaleN = 1.0;
+ opt->scaleM = 1.0; opt->scaleN = 1.0; opt->zeroN = 0; opt->zeroV = 1.0;
ierr = PetscOptionsBegin(comm, "", "Stokes problem options", "");CHKERRQ(ierr);
ierr = PetscOptionsIntArray("-dim", "list of dimension extent", "stokes.C", ctx->dim, &ctx->numDims, &flag);CHKERRQ(ierr);
if (!flag) { ctx->numDims = 2; ctx->dim[0] = 8; ctx->dim[1] = 6; }
@@ -409,6 +409,8 @@ PetscErrorCode StokesProcessOptions(StokesCtx *ctx)
ierr = PetscOptionsInt("-cont", "number of continuations", "stokes.C", opt->cont, &opt->cont, PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-scaleM", "scaling factor for Mixed condition", "stokes.C", opt->scaleM, &opt->scaleM, PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsReal("-scaleN", "scaling factor for Neumann condition", "stokes.C", opt->scaleN, &opt->scaleN, PETSC_NULL);CHKERRQ(ierr);
+ ierr = PetscOptionsInt("-zeroN", "number of rows to zero", "stokes.C", opt->zeroN, &opt->zeroN, PETSC_NULL);CHKERRQ(ierr);
+ ierr = PetscOptionsReal("-zeroV", "value to put on diagonal", "stokes.C", opt->zeroV, &opt->zeroV, PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsEnd();
ierr = PetscPrintf(comm, "Stokes problem");CHKERRQ(ierr);
@@ -1401,6 +1403,20 @@ PetscErrorCode StokesPCSetUp1(void *void_ctx)
lump[row[i]] += M[i][j];
}
}
+
+ if (ctx->options->zeroN) { // filter the rows and columns to produce a symmetric matrix with some points fixed
+ for (int i=0; i<d*N; i++) {
+ if (row[i] < ctx->options->zeroN) {
+ for (int j=0; j<d*N; j++) {
+ if (col[j] == row[i]) {
+ A[i][j] = ctx->options->zeroV;
+ } else {
+ A[i][j] = A[j][i] = 0.0;
+ }
+ }
+ }
+ }
+ }
ierr = MatSetValues(ctx->MatVVPC, d*N, row, d*N, col, &A[0][0], ADD_VALUES);CHKERRQ(ierr);
if (ctx->options->debug > 1) {
ierr = MatAssemblyBegin(ctx->MatVVPC, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

0 comments on commit c4ec55a

Please sign in to comment.