Permalink
Browse files

Poisson solver.

darcs-hash:20080303223426-40437-d68e6bd220058e767326782d29634ec83d089f16.gz
  • Loading branch information...
1 parent b2d95fb commit 92d68dd3b8a7929e26dd6f77da703d9554b7b838 @jedbrown committed Mar 3, 2008
Showing with 358 additions and 24 deletions.
  1. +4 −1 Makefile
  2. +26 −9 cheb.c
  3. +15 −14 chebyshev.c
  4. +313 −0 poisson.c
View
@@ -1,7 +1,7 @@
ALL : ex15 shell
CFLAGS= -std=c99
-include ${PETSC_DIR}/bmake/common/base
+include ${PETSC_DIR}/conf/base
nk : nk.o chkopts
${CLINKER} -o nk nk.o ${PETSC_LIB}
@@ -14,3 +14,6 @@ shell : shell.o chkopts
cheb : cheb.o chebyshev.o chkopts
${CLINKER} -o cheb cheb.o chebyshev.o -lfftw3 ${PETSC_LIB}
+
+poisson : poisson.o chebyshev.o chkopts
+ ${CLINKER} -o poisson poisson.o chebyshev.o -lfftw3 ${PETSC_LIB}
View
35 cheb.c
@@ -17,7 +17,7 @@ int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A;
- Vec x2,b2,u2; /* approx solution, RHS, exact solution */
+ Vec x2,b2,u2,du2; /* approx solution, RHS, exact solution */
Mat A2;
KSP ksp; /* linear solver context */
PC pc; /* preconditioner context */
@@ -45,16 +45,17 @@ int main(int argc,char **args)
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
ierr = MatCreateChebD1(PETSC_COMM_WORLD, x, b,
- FFTW_ESTIMATE | FFTW_PRESERVE_INPUT, &A); CHKERRQ(ierr);
+ FFTW_ESTIMATE, &A); CHKERRQ(ierr);
int dims[] = { m, n, p };
// ierr = VecCreate(PETSC_COMM_WORLD, &u2);CHKERRQ(ierr);
// ierr = VecSetSizes(u2, PETSC_DECIDE, m * n);CHKERRQ(ierr);
// ierr = VecSetFromOptions(u2);CHKERRQ(ierr);
ierr = VecCreateSeq(PETSC_COMM_WORLD, m * n * p, &u2); CHKERRQ(ierr);
ierr = VecDuplicate(u2, &b2);CHKERRQ(ierr);
- ierr = VecDuplicate(b2, &x2);CHKERRQ(ierr);
- ierr = MatCreateCheb(PETSC_COMM_WORLD, 3, d, dims, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT,
+ ierr = VecDuplicate(u2, &du2);CHKERRQ(ierr);
+ ierr = VecDuplicate(u2, &x2);CHKERRQ(ierr);
+ ierr = MatCreateCheb(PETSC_COMM_WORLD, 3, d, dims, FFTW_ESTIMATE,
x2, b2, &A2); CHKERRQ(ierr);
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
@@ -70,21 +71,28 @@ int main(int argc,char **args)
ierr = VecRestoreArray(u, &a); CHKERRQ(ierr);
// 2-D solution function
+ double *e;
ierr = VecGetArray(u2, &a); CHKERRQ(ierr);
+ ierr = VecGetArray(du2, &e); CHKERRQ(ierr);
for (int i=0; i < m; i++) {
double x = (m==1) ? 0 : cos (i * PI / (m-1));
for (int j=0; j < n; j++) {
double y = (n==1) ? 0 : cos (j * PI / (n-1));
for (int k=0; k < p; k++) {
double z = (p==1) ? 0 : cos (k * PI / (p-1));
a[(i*n + j) * p + k] = exp(x) + exp(y) + exp(z);
+ switch (d) {
+ case 0: e[(i*n + j) * p + k] = exp(x); break;
+ case 1: e[(i*n + j) * p + k] = exp(y); break;
+ case 2: e[(i*n + j) * p + k] = exp(z); break;
+ }
}
}
}
ierr = VecRestoreArray(u2, &a); CHKERRQ(ierr);
+ ierr = VecRestoreArray(du2, &e); CHKERRQ(ierr);
ierr = MatMult(A, u, b); CHKERRQ(ierr);
- ierr = MatMult(A2, u2, b2); CHKERRQ(ierr);
// ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ierr = VecCopy(b, x); CHKERRQ(ierr);
@@ -94,11 +102,20 @@ int main(int argc,char **args)
//ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A iterations %D\n",norm,its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A\n",norm);CHKERRQ(ierr);
+ ierr = MatMult(A2, u2, b2); CHKERRQ(ierr);
+ ierr = VecCopy(b2, x2); CHKERRQ(ierr);
+
+ ierr = VecAXPY(x2,none,du2);CHKERRQ(ierr);
+ ierr = VecNorm(x2,NORM_INFINITY,&norm);CHKERRQ(ierr);
+ // ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
+ //ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A iterations %D\n",norm,its);CHKERRQ(ierr);
+ ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %A\n",norm);CHKERRQ(ierr);
+
/* /\* ierr = VecView(x, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr); *\/ */
- printf("\n");
- ierr = VecView(b, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
- printf("\n");
- ierr = VecView(b2, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr);
+ /* printf("\n"); */
+ /* ierr = VecView(b, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr); */
+ /* printf("\n"); */
+ /* ierr = VecView(b2, PETSC_VIEWER_STDOUT_SELF); CHKERRQ(ierr); */
/*
Free work space. All PETSc objects should be destroyed when they
View
@@ -123,9 +123,9 @@ PetscErrorCode MatCreateCheb(MPI_Comm comm, int rank, int tr, int *dim, unsigned
ierr = VecGetArray(vx, &x); CHKERRQ(ierr);
ierr = VecGetArray(vy, &y); CHKERRQ(ierr);
const fftw_r2r_kind redft00 = FFTW_REDFT00, rodft00 = FFTW_RODFT00;
- c->p_forward = fftw_plan_guru_r2r(1, &(c->tdim), c->rank-1, c->dim, x, c->work, &redft00, flag);
+ c->p_forward = fftw_plan_guru_r2r(1, &(c->tdim), c->rank-1, c->dim, x, c->work, &redft00, flag | FFTW_PRESERVE_INPUT);
fftw_iodim tdim1 = c->tdim; tdim1.n -= 2; // we need a modified tdim to come back.
- c->p_backward = fftw_plan_guru_r2r(1, &tdim1, c->rank-1, c->dim, c->work + tdim1.os, y + tdim1.is, &rodft00, flag);
+ c->p_backward = fftw_plan_guru_r2r(1, &tdim1, c->rank-1, c->dim, c->work + tdim1.os, y + tdim1.is, &rodft00, flag | FFTW_DESTROY_INPUT);
ierr = VecRestoreArray(vx, &x); CHKERRQ(ierr);
ierr = VecRestoreArray(vy, &y); CHKERRQ(ierr);
@@ -153,38 +153,39 @@ PetscErrorCode ChebMult(Mat A, Vec vx, Vec vy) {
fftw_execute_r2r(c->p_forward, x, c->work);
+ double N = (double)n;
ierr = PetscMemzero(ind, (c->rank - 1) * sizeof(int)); CHKERRQ(ierr);
int offset = 0;
for (bool done = false; !done; ) {
+ int ix0 = offset;
+ int ixn = offset + n * c->tdim.is;
+ y[ix0] = 0.0;
+ y[ixn] = 0.0;
+ double s = 1.0;
for (int i = 1; i < n; i++) {
int ix = offset + i * c->tdim.is;
- c->work[ix] *= (double)i;
+ double I = (double)i;
+ c->work[ix] *= I;
+ y[ix0] += I * c->work[ix];
+ y[ixn] += s * I * c->work[ix];
+ s = -s;
}
+ y[ix0] = 0.5 * c->work[ixn] * N + y[ix0] / n;
+ y[ixn] = y[ixn] / N + 0.5 * s * N * c->work[ixn];
perform_carry(c, ind, &offset, &done);
}
fftw_execute_r2r(c->p_backward, c->work + c->tdim.os, y + c->tdim.is);
- double N = (double)n;
double pin = PI / N;
ierr = PetscMemzero(ind, (c->rank - 1) * sizeof(int)); CHKERRQ(ierr);
offset = 0;
for (bool done = false; !done; ) {
- int ix0 = offset;
- int ixn = offset + n * c->tdim.is;
- y[ix0] = 0.0;
- y[ixn] = 0.0;
- double s = 1.0;
for (int i = 1; i < n; i++) {
int ix = offset + i * c->tdim.is;
double I = (double)i;
y[ix] /= 2 * n * sqrt(1.0 - PetscSqr(cos(I * pin)));
- y[ix0] += I * c->work[ix];
- y[ixn] += s * I * c->work[ix];
- s = -s;
}
- y[ix0] = 0.5 * c->work[ixn] * N + y[ix0] / n;
- y[ixn] = y[ixn] / N + 0.5 * s * N * c->work[ixn];
perform_carry(c, ind, &offset, &done);
}
Oops, something went wrong.

0 comments on commit 92d68dd

Please sign in to comment.