-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1
Dear PETSc users,
I'm trying to write a Poisson solver in Fortran using the
KSPSetComputeOperators etc. framework. When debugging this, I ended up
modifying ksp/ex22.f so that it matches ksp/ex34.c. The difference
between these is that ex34.c has a non-constant RHS and Neumann BCs,
which is closer to what I want.
Now, when I run these two programs I get the following:
./ex22f_mod
Residual, L2 norm 1.007312
Error, sup norm 0.020941
Error, L1 norm 10.687882
Error, L2 norm 0.340425
./ex34
Residual, L2 norm 1.07124e-05
Error, sup norm 0.0209405
Error, L1 norm 0.00618512
Error, L2 norm 0.000197005
but when I MatView/VecView the matrix and RHS for each example (write
them to file), there is no difference. I have attached ex22f_mod.F90,
any suggestions on what the error is?
Once this example works, you can of course include it with PETSc.
Regards,
Åsmund
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.22 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/
iQEcBAEBAgAGBQJTPWWZAAoJED+FDAHgGz19aNYIANyq9k3a5kHzTlAdybkZCYDw
yKtDE5l/4iWYmNL49FH8AocHVcRivLaeJG5CGKqySFtZUXOlC9DM7rn4UmuQecni
dgIQuTk0Ym+OJccHyT5xxnebpFVNrIOTpInfQaDW6dTyeL1svAMeHqslKaGepySL
q/cODbgNDYgl6uumB+POMZevtlM6HPhl/1m7HofcHC9upvTRjSPqP1cg+kg+/8m2
Fdie/X7PCBfShrAys94kNXNcwtbO7taauphkQGMfyl0gUd+lFATG6zrEZdDqSFlV
c44GFFbxW/SRIDBXdOeX9/cy75KW5do1Sildwb6R4H/i7t6/hCJUJuss7FHjmLc=
=tV3N
-----END PGP SIGNATURE-----
!
! Laplacian in 3D. Modeled by the partial differential equation
!
! Laplacian u = 1,0 < x,y,z < 1,
!
! with boundary conditions
!
! u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
!
! This uses multigrid to solve the linear system
program main
implicit none
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Include files
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! petscsys.h - base PETSc routines petscvec.h - vectors
! petscmat.h - matrices
! petscksp.h - Krylov subspace methods petscpc.h - preconditioners
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscmat.h>
#include <finclude/petscpc.h>
#include <finclude/petscksp.h>
#include <finclude/petscdmda.h>
#include <finclude/petscdmda.h90>
PetscErrorCode ierr
DM da
KSP ksp
Vec x,b,r
PetscInt i1,i3
PetscInt ctx
PetscScalar r1, norm,Hx,Hy,hz
MatNullSpace nullspace
character(len=255) :: Fstr
PetscScalar, pointer :: array(:,:,:)
PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs
interface
subroutine ComputeRHS(ksp,rhs,ctx,ierr)
KSP, intent(in) :: ksp
Vec, intent(inout) :: rhs
PetscInt, intent(in) :: ctx
PetscErrorCode, intent(inout) :: ierr
end subroutine ComputeRHS
!
subroutine ComputeMatrix(ksp,JJ,A,str,ctx,ierr)
KSP, intent(in) :: ksp
Mat, intent(in) :: A,JJ
MatStructure, intent(in) :: str
PetscInt, intent(in) :: ctx
PetscErrorCode, intent(in) :: ierr
end subroutine ComputeMatrix
end interface
call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
i3 = -12
i1 = 1
call KSPCreate(MPI_COMM_WORLD,ksp,ierr)
call DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE, &
& DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE, &
& DMDA_STENCIL_STAR,i3,i3,i3,PETSC_DECIDE,PETSC_DECIDE, &
& PETSC_DECIDE,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,da,ierr)
call DMDASetInterpolationType(da, DMDA_Q0, ierr)
call KSPSetDM(ksp,da,ierr)
! call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,nullspace,ierr)
! call KSPSetNullSpace(ksp,nullspace,ierr)
call KSPSetComputeRHS(ksp,ComputeRHS,ctx,ierr)
call KSPSetComputeOperators(ksp,ComputeMatrix,ctx,ierr)
call KSPSetFromOptions(ksp,ierr)
call KSPSolve(ksp,PETSC_NULL_OBJECT,PETSC_NULL_OBJECT,ierr)
call KSPGetSolution(ksp,x,ierr)
call KSPGetRhs(ksp,b,ierr)
call KSPGetOperators(ksp,PETSC_NULL_OBJECT,J,PETSC_NULL_OBJECT,ierr)
call VecDuplicate(b,r,ierr)
call MatMult(J,x,r,ierr)
call VecAXPY(r,-1.0,b,ierr)
call VecNorm(r,NORM_2,norm,ierr)
write(Fstr,101) norm
call PetscPrintf(PETSC_COMM_WORLD,Fstr,ierr)
call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,ierr)
Hx = 1.d0 / (mx)
Hy = 1.d0 / (my)
Hz = 1.d0 / (mz)
call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
call DMDAVecGetArrayF90(da,x,array,ierr)
r1=1.d0
do k=zs,zs+zm-1
do j=ys,ys+ym-1
do i=xs,xs+xm-1
array(i,j,k) = array(i,j,k) - &
cos(2*PETSC_PI*((r1*i+0.5)*Hx))* &
cos(2*PETSC_PI*((r1*j+0.5)*Hy))* &
cos(2*PETSC_PI*((r1*k+0.5)*Hz))
end do
end do
end do
call DMDAVecRestoreArrayF90(da,x,array,ierr)
call VecAssemblyBegin(x,ierr)
call VecAssemblyEnd(x,ierr)
call VecNorm(x,NORM_INFINITY,norm,ierr)
write(Fstr,102) norm
call PetscPrintf(PETSC_COMM_WORLD,Fstr,ierr)
call VecNorm(x,NORM_1,norm,ierr)
write(Fstr,103) norm
call PetscPrintf(PETSC_COMM_WORLD,Fstr,ierr)
call VecNorm(x,NORM_2,norm,ierr)
write(Fstr,104) norm
call PetscPrintf(PETSC_COMM_WORLD,Fstr,ierr)
call KSPDestroy(ksp,ierr)
call DMDestroy(da,ierr)
call PetscFinalize(ierr)
101 format(1x,"Residual, L2 norm ",F16.6,"\n")
102 format(1x,"Error, sup norm ",F16.6,"\n")
103 format(1x,"Error, L1 norm ",F16.6,"\n")
104 format(1x,"Error, L2 norm ",F16.6,"\n")
end
subroutine ComputeRHS(ksp,b,ctx,ierr)
implicit none
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscdmda.h>
#include <finclude/petscdmda.h90>
KSP ksp
Vec b
PetscInt ctx
PetscErrorCode ierr
PetscInt mx,my,mz
PetscScalar Hx,Hy,Hz,r1
DM da
PetscInt i,j,k,xm,ym,zm,xs,ys,zs
PetscScalar, pointer :: array(:,:,:)
MatNullSpace nullspace
PetscViewer viewer
call KSPGetDM(ksp,da,ierr)
call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& PETSC_NULL_INTEGER,ierr)
Hx=1.d0/(mx)
Hy=1.d0/(my)
Hz=1.d0/(mz)
call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr)
call DMDAVecGetArrayF90(da,b,array,ierr)
r1=1.d0
do k=zs,zs+zm-1
do j=ys,ys+ym-1
do i=xs,xs+xm-1
array(i,j,k) = 12*PETSC_PI*PETSC_PI* &
cos(2*PETSC_PI*((r1*i+0.5)*Hx))* &
cos(2*PETSC_PI*((r1*j+0.5)*Hy))* &
cos(2*PETSC_PI*((r1*k+0.5)*Hz)) * Hx*Hy*Hz
end do
end do
end do
call DMDAVecRestoreArrayF90(da,b,array,ierr)
call VecAssemblyBegin(b,ierr)
call VecAssemblyEnd(b,ierr)
call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,nullspace,ierr)
call MatNullSpaceRemove(nullspace,b,PETSC_NULL_OBJECT,ierr)
call MatNullSpaceDestroy(nullspace,ierr)
call PetscViewerASCIIOpen(PETSC_COMM_WORLD,"vec22.out",viewer,ierr)
call VecView(b,viewer,ierr)
return
end
!---------------------------------------------------------------------
subroutine ComputeMatrix(ksp,JJ,A,str,ctx,ierr)
implicit none
#include <finclude/petscsys.h>
#include <finclude/petscvec.h>
#include <finclude/petscviewer.h>
#include <finclude/petscdmda.h>
#include <finclude/petscdmda.h90>
#include "finclude/petscmat.h"
#include "finclude/petscksp.h"
#include "finclude/petscpc.h"
!
KSP, intent(in) :: ksp
Mat, intent(in) :: A,JJ
MatStructure, intent(in) :: str
PetscInt, intent(in) :: ctx
PetscErrorCode, intent(in) :: ierr
!
DM da
PetscInt i,j,k,mx,my,mz,im
PetscInt jm,km,is,js,ks,i1,i7,MSi,MSj,MSk
PetscInt :: num,numi,numj,numk
PetscScalar v(7)
PetscScalar HxHydHz,HyHzdHx,HxHzdHy,dx,dy,dz
MatStencil row(4),col(4,7)
MatNullSpace nullspace
PetscViewer viewer
!
i1 = 1
i7 = 7
MSi=MatStencil_i
MSj=MatStencil_j
MSk=MatStencil_k
!
call KSPGetDM(ksp,da,ierr)
call DMDAGetCorners(da,is,js,ks,im,jm,km,ierr)
im=is+im-1
jm=js+jm-1
km=ks+km-1
call DMDAGetInfo(da,PETSC_NULL_INTEGER,mx,my,mz,&
PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,&
PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,&
PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)
!
dx=1.d0/(mx)
dy=1.d0/(my)
dz=1.d0/(mz)
HxHydHz = dx*dy/dz
HxHzdHy = dx*dz/dy
HyHzdHx = dy*dz/dx
!
! Set the matrix values
do k=ks,km
do j=js,jm
do i=is,im
row(MSi) = i
row(MSj) = j
row(MSk) = k
if (i==0 .or. j==0 .or. k==0 .or. &
i==mx-1 .or. j==my-1 .or. k==mz-1) then
! BC's
num=1
numi=0
numj=0
numk=0
! Lots of if's to check the edges.
!
! Eg. if we are not at the bottom, include the bottom stencil part
if (k/=0) then
v(num) = -HxHydHz
col(MSi,num) = i
col(MSj,num) = j
col(MSk,num) = k-1
num=num+1
numk=numk+1
end if
if (j/=0) then
v(num) = -HxHzdHy
col(MSi,num) = i
col(MSj,num) = j-1
col(MSk,num) = k
num=num+1
numj=numj+1
end if
if (i/=0) then
v(num) = -HyHzdHx
col(MSi,num) = i-1
col(MSj,num) = j
col(MSk,num) = k
num=num+1
numi=numi+1
end if
if (i/=mx-1) then
v(num) = -HyHzdHx
col(MSi,num) = i+1
col(MSj,num) = j
col(MSk,num) = k
num=num+1
numi=numi+1
end if
if (j/=my-1) then
v(num) = -HxHzdHy
col(MSi,num) = i
col(MSj,num) = j+1
col(MSk,num) = k
num=num+1
numj=numj+1
end if
if (k/=mz-1) then
v(num) = -HxHydHz
col(MSi,num) = i
col(MSj,num) = j
col(MSk,num) = k+1
num=num+1
numk=numk+1
end if
! Central point
v(num) = numk*HxHydHz + numj*HxHzdHy + numi*HyHzdHx
col(MSi,num) = i
col(MSj,num) = j
col(MSk,num) = k
call MatSetValuesStencil(A,1,row,num,col,v,INSERT_VALUES,ierr)
CHKERRQ(ierr)
else
! Standard stencil away from boundaries
col(MSi,:)=i
col(MSj,:)=j
col(MSk,:)=k
v(1) = -HxHydHz; col(MSk,1) = k-1
v(2) = -HxHzdHy; col(MSj,2) = j-1
v(3) = -HyHzdHx; col(MSi,3) = i-1
v(4) = 2.0*(HyHzdHx + HxHzdHy + HxHydHz)
v(5) = -HyHzdHx; col(MSi,5) = i+1
v(6) = -HxHzdHy; col(MSj,6) = j+1
v(7) = -HxHydHz; col(MSk,7) = k+1
! Set these values into the matrix
call MatSetValuesStencil(A,i1,row,i7,col,v,INSERT_VALUES,ierr)
CHKERRQ(ierr)
end if
end do
end do
end do
!
! Assemble the A matrix and right-hand-side vector
call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
CHKERRQ(ierr)
call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
CHKERRQ(ierr)
call MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,nullspace,ierr)
call MatSetNullSpace(A,nullspace,ierr)
call MatNullSpaceDestroy(nullspace,ierr)
call PetscViewerASCIIOpen(PETSC_COMM_WORLD,"mat22.out",viewer,ierr)
call MatView(A,viewer,ierr)
end subroutine ComputeMatrix
