[petsc-users] An issue of order of accuracy.
Alan Z. Wei
zhenglun.wei at gmail.com
Wed Dec 4 17:25:12 CST 2013
Dear all,
I hope you had a great Thanksgiving.
Currently, I tested the order of accuracy for
/src/ksp/ksp/tutorial/example/ex45.c. Since the 2nd-order discretization
is used in this program and ksp solver is converged to 10^-7, I expected
that the solution should provides a 2nd-order in L2 norm. However, as I
tested (even with a Laplace equation), the L2 norm slope is much less
than 2. Sometime, if the grid size is reduced, the L2 norm increases.
Could anyone help me about this issue, please?
Here is the L2 norm outputted:
Grid L2 norm (10^-8)
0.05 4.36242
0.025 2.20794
0.0125 7.02749
0.00625 12.64
Once the grid size is reduced to half, the number of the grid will be
multiplied by 8 in order to keep the same size of the computational domain.
The code is also attached. It is from ex45.c with very little
modifications.
thanks in advance,
Alan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131204/9de2c2e2/attachment.html>
-------------- next part --------------
static char help[] = "Solves 3D Laplacian using multigrid.\n\n";
#include <petscksp.h>
#include <petscdmda.h>
const double dx = 0.0125;
const double dy = 0.0125;
const double dz = 0.0125;
const int imax = 40;
const int jmax = 40;
const int kmax = 40;
extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
extern PetscErrorCode ComputeInitialGuess(KSP,Vec,void*);
extern PetscErrorCode PostP(DM, Vec, Vec);
PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
PetscErrorCode ierr;
KSP ksp;
PetscReal norm;
DM da;
Vec x,b,r;
Mat A;
double t1, t2;
PetscInitialize(&argc,&argv,(char*)0,help);
t1 = MPI_Wtime();
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-imax,-jmax,-kmax,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,da);CHKERRQ(ierr);
ierr = KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,NULL);CHKERRQ(ierr);
ierr = KSPSetComputeRHS(ksp,ComputeRHS,NULL);CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp,ComputeMatrix,NULL);CHKERRQ(ierr);
ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,NULL,NULL);CHKERRQ(ierr);
ierr = KSPGetSolution(ksp,&x);CHKERRQ(ierr);
ierr = KSPGetRhs(ksp,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&r);CHKERRQ(ierr);
ierr = KSPGetOperators(ksp,&A,NULL,NULL);CHKERRQ(ierr);
ierr = MatMult(A,x,r);CHKERRQ(ierr);
ierr = VecAXPY(r,-1.0,b);CHKERRQ(ierr);
ierr = VecNorm(r,NORM_2,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);CHKERRQ(ierr);
ierr = KSPGetDM(ksp, &da);CHKERRQ(ierr);
ierr = PostP(da, x, b);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
t2 = MPI_Wtime();
printf("Total Time Elapsed: %f\n", t2-t1);
ierr = PetscFinalize();
return 0;
}
#undef __FUNCT__
#define __FUNCT__ "ComputeRHS"
PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
{
PetscErrorCode ierr;
PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs, mm,nn,pp;
DM dm;
PetscScalar Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
PetscScalar ***barray;
PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&dm);CHKERRQ(ierr);
ierr = DMDAGetInfo(dm,0,&mx,&my,&mz,&mm,&nn,&pp,0,0,0,0,0,0);CHKERRQ(ierr);
printf("mx = %d, my = %d, mz =%d, mm = %d, nn = %d, pp = %d\n", mx, my, mz, mm, nn, pp);
// Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
Hx = dx; Hy = dy; Hz = dz;
HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
ierr = DMDAVecGetArray(dm,b,&barray);CHKERRQ(ierr);
for (k=zs; k<zs+zm; k++) {
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
barray[k][j][i] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
} else {
barray[k][j][i] = 0.;
}
}
}
}
ierr = DMDAVecRestoreArray(dm,b,&barray);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "ComputeInitialGuess"
PetscErrorCode ComputeInitialGuess(KSP ksp,Vec b,void *ctx)
{
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = VecSet(b,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "ComputeMatrix"
PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,MatStructure *stflg,void *ctx)
{
DM da;
PetscErrorCode ierr;
PetscScalar v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
MatStencil row,col[7];
PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
ierr = DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
// Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
Hx = dx; Hy = dy; Hz = dz;
HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr);
for (k=zs; k<zs+zm; k++) {
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.i = i; row.j = j; row.k = k;
if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
v[0] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
ierr = MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr);
} else {
v[0] = -HxHydHz;col[0].i = i; col[0].j = j; col[0].k = k-1;
v[1] = -HxHzdHy;col[1].i = i; col[1].j = j-1; col[1].k = k;
v[2] = -HyHzdHx;col[2].i = i-1; col[2].j = j; col[2].k = k;
v[3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);col[3].i = row.i; col[3].j = row.j; col[3].k = row.k;
v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
ierr = MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr);
}
}
}
}
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
*stflg = SAME_NONZERO_PATTERN;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "PostP"
PetscErrorCode PostP(DM da, Vec x, Vec b) {
PetscErrorCode ierr;
PetscScalar ***DpArray;
double L2Norm;
FILE *fpr;
PetscFunctionBeginUser;
fpr = fopen("Check.dat", "w");
L2Norm = 0.;
ierr = DMDAVecGetArray(da, x, &DpArray);CHKERRQ(ierr);
for (k=zs+1; k<zs+zm-1; k++) {
for (j=ys+1; j<ys+ym-1; j++) {
for (i=xs+1; i<xs+xm-1; i++) {
// L2Norm += pow(DpArray[k][j][i] - (1000.*(i*dx)*(i*dx)/2.+(i*dx)), 2.0);
L2Norm += pow(DpArray[k][j][i] - 1., 2.0);
fprintf(fpr, "%f %f %f %.10f \n", i*dx, j*dy, k*dz, DpArray[k][j][i]);
}
}
}
fclose(fpr);
L2Norm = pow(L2Norm/(mx-2)/(my-2)/(mz-2), 0.5);
printf("L2Norm = %e\n", L2Norm);
ierr = DMDAVecRestoreArray(da, x, &DpArray);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
-------------- next part --------------
rm -f Check.dat
make ex45
#/bio/work1/zlwei/PETSc/petsc-dev/linux-gnu-c-nodebug/bin/mpiexec -np 128 ./ex45 -pc_type mg -ksp_type cg -da_refine 7 -pc_mg_galerkin -ksp_rtol 1.0e-7 -mg_levels_pc_type jacobi -mg_levels_ksp_type chebyshev
#/bio/work1/zlwei/PETSc/petsc-dev/linux-gnu-c-nodebug/bin/mpiexec -np 32 ./ex45 -pc_type mg -ksp_type cg -da_refine 5 -pc_mg_galerkin -ksp_rtol 1.0e-7 -mg_levels_pc_type jacobi -mg_levels_ksp_type chebyshev -mg_coarse_ksp_type cg -mg_coarse_pc_type jacobi > out &
#/bio/work1/zlwei/PETSc/petsc-dev/linux-gnu-c-nodebug/bin/mpiexec -np 32 ./ex45 -pc_type mg -ksp_type cg -da_refine 5 -pc_mg_galerkin -ksp_rtol 1.0e-7 -mg_levels_pc_type jacobi -mg_levels_ksp_type chebyshev -dm_view -log_summary -pc_mg_log -pc_mg_monitor -ksp_view -ksp_monitor -pc_mg_log> out &
/bio/work1/zlwei/PETSc/petsc-dev/linux-gnu-c-nodebug/bin/mpiexec -np 1 ./ex45 -pc_type gamg -ksp_type gmres -pc_gamg_agg_nsmooths 1 -mg_levels_ksp_max_it 1 -mg_levels_ksp_type richardson -ksp_rtol 1.0e-7
More information about the petsc-users
mailing list