! Concepts: KSP^solving a system of linear equations ! Concepts: KSP^Laplacian, 3d ! Processors: n ! !Laplacian in 3D. Modeled by the partial differential equation ! ! div grad u = f, 0 < x,y,z < 1, ! !with pure Neumann boundary conditions ! ! u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. ! !The functions are cell-centered ! !This uses multigrid to solve the linear system ! ! Contributed by Jianming Yang ! !changed by peterspy to fortran program main #include #include #include use petscdm use petscdmda use petscksp implicit none DM da KSP ksp Vec x,b,r Mat JJ PetscErrorCode ierr PetscReal norm,neg_one PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,d,ctx PetscScalar Hx,Hy,Hz,array(12,12,12) external ComputeRHS,ComputeMatrix neg_one=-1.0 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) if (ierr .ne. 0) then print*,'Unable to initialize PETSc' stop endif call KSPCreate(MPI_COMM_WORLD,ksp,ierr) call DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, & DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,12,12,12, & PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr) call DMSetFromOptions(da,ierr) call DMSetUp(da,ierr) call DMDASetInterpolationType(da,DMDA_Q0,ierr) call KSPSetDM(ksp,da,ierr) call KSPSetComputeRHS(ksp,ComputeRHS,ctx,ierr) call KSPSetComputeOperators(ksp,ComputeMatrix,ctx,ierr) call KSPSetFromOptions(ksp,ierr) call KSPSolve(ksp,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr) call KSPGetSolution(ksp,x,ierr) call KSPGetRhs(ksp,b,ierr) call KSPGetOperators(ksp,PETSC_NULL_MAT,JJ,ierr) call VecDuplicate(b,r,ierr) call MatMult(JJ,x,r,ierr) call VecAXPY(r,neg_one,b,ierr) call VecNorm(r,NORM_2,norm,ierr) call PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",real(norm),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.0/real(mx-1) Hy=1.0/real(my-1) Hz=1.0/real(mz-1) call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr) do 10 k=zs,zs+zm-1 do 10 j=ys,ys+ym-1 do 10 i=xs,xs+xm-1 array(i,j,k)=cos(2*PETSC_PI*((real(i)+0.5)*Hx))* & cos(2*PETSC_PI*((real(j)+0.5)*Hy))* & cos(2*PETSC_PI*((real(k)+0.5)*Hz)) 10 continue call VecAssemblyBegin(x,ierr) call VecAssemblyEnd(x,ierr) call VecNorm(x,NORM_INFINITY,norm,ierr) call PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",real(norm),ierr) call VecNorm(x,NORM_1,norm,ierr) call PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",real(norm/(real(mx)*real(my)*real(mz))),ierr) call VecNorm(x,NORM_2,norm,ierr) call PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",real(norm/(real(mx)*real(my)*real(mz))),ierr) call VecDestroy(r,ierr) call KSPDestroy(ksp,ierr) call DMDestroy(da,ierr) call PetscFinalize(ierr) endprogram main subroutine ComputeRHS(ksp,b,ctx,ierr) use petscksp implicit none PetscErrorCode ierr PetscInt d,i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,ctx PetscScalar Hx,Hy,Hz,array(12,12,12) Vec b KSP ksp DM da 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.0/real(mx) Hy=1.0/real(my) Hz=1.0/real(mz) call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr) do 20 k=zs,zs+zm-1 do 20 j=ys,ys+ym-1 do 20 i=xs,xs+xm-1 array(i,j,k)=12*PETSC_PI*PETSC_PI*Hx*Hy*Hz* & cos(2*PETSC_PI*((real(i)+0.5)*Hx))* & cos(2*PETSC_PI*((real(j)+0.5)*Hy))* & cos(2*PETSC_PI*((real(k)+0.5)*Hz)) 20 continue call VecAssemblyBegin(b,ierr) call VecAssemblyEnd(b,ierr) ! /* force right hand side to be consistent for singular matrix */ ! /* note this is really a hack, normally the model would provide you with a !consistent right handside */ return endsubroutine ComputeRHS subroutine ComputeMatrix(ksp,JJ,jac,ctx,ierr) use petscksp implicit none PetscErrorCode ierr Mat jac,JJ KSP ksp PetscInt i,j,k,d,mx,my,mz,xm,ym,zm,xs,ys,zs,num,numi,numj,numk,ctx PetscScalar v(7),Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz MatStencil row(4),col(4,7) DM da 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.0/real(mx) Hy=1.0/real(my) Hz=1.0/real(mz) HyHzdHx=Hy*Hz/Hx HxHzdHy=Hx*Hz/Hy HxHydHz=Hx*Hy/Hz call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr) do 30 k=zs,zs+zm-1 do 30 j=ys,ys+ym-1 do 30 i=xs,xs+xm-1 row(MatStencil_i) = i row(MatStencil_j) = j row(MatStencil_k) = k if (i.eq.0 .or. j.eq.0 .or. k.eq.0 .or. i.eq.mx-1 .or. j.eq.my-1 .or. k.eq.mz-1) then num=1;numi=0;numj=0;numk=0 if (k.ne.0) then v(num) = -HxHydHz col(MatStencil_i,num) = i col(MatStencil_j,num) = j col(MatStencil_k,num) = k-1 num=num+1;numk=numk+1 endif if (j.ne.0) then v(num) = -HxHzdHy col(MatStencil_i,num) = i col(MatStencil_j,num) = j-1 col(MatStencil_k,num) = k num=num+1;numj=numj+1 endif if (i.ne.0) then v(num) = -HyHzdHx col(MatStencil_i,num) = i-1 col(MatStencil_j,num) = j col(MatStencil_k,num) = k num=num+1;numi=numi+1 endif if (i.ne.mx-1) then v(num) = -HyHzdHx col(MatStencil_i,num) = i+1 col(MatStencil_j,num) = j col(MatStencil_k,num) = k num=num+1;numi=numi+1 endif if (j.ne.my-1) then v(num) = -HxHzdHy col(MatStencil_i,num) = i col(MatStencil_j,num) = j+1 col(MatStencil_k,num) = k num=num+1;numj=numj+1 endif if (k.ne.mz-1) then v(num) = -HxHydHz col(MatStencil_i,num) = i col(MatStencil_j,num) = j col(MatStencil_k,num) = k+1 num=num+1;numk=numk+1 endif v(num)=real(numk)*HxHydHz+real(numj)*HxHzdHy+real(numi)*HyHzdHx col(MatStencil_i,num) = i col(MatStencil_j,num) = j col(MatStencil_k,num) = k num=num+1 call MatSetValuesStencil(jac,1,row,num,col,v,INSERT_VALUES,ierr) else v(1)=-HxHydHz; col(MatStencil_i,1)=i; col(MatStencil_j,1)=j; col(MatStencil_k,1)=k-1; v(2)=-HxHzdHy; col(MatStencil_i,2)=i; col(MatStencil_j,2)=j-1; col(MatStencil_k,2)=k; v(3)=-HyHzdHx; col(MatStencil_i,3)=i-1; col(MatStencil_j,3)=j; col(MatStencil_k,3)=k; v(4)=2.0*(HyHzdHx + HxHzdHy + HxHydHz); col(MatStencil_i,4)=i; col(MatStencil_j,4)=j; col(MatStencil_k,4)=k; v(5)=-HyHzdHx; col(MatStencil_i,5)=i+1; col(MatStencil_j,5)=j; col(MatStencil_k,5)=k; v(6)=-HxHzdHy; col(MatStencil_i,6)=i; col(MatStencil_j,6)=j+1; col(MatStencil_k,6)=k; v(7)=-HxHydHz; col(MatStencil_i,7)=i; col(MatStencil_j,7)=j; col(MatStencil_k,7)=k+1; call MatSetValuesStencil(jac,1,row,7,col,v,INSERT_VALUES,ierr) endif 30 continue call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr) return endsubroutine ComputeMatrix !/*TEST ! ! build: ! requires: !complex !single ! ! test: ! args: -pc_type mg -pc_mg_type full -ksp_type fgmres -ksp_monitor_short -pc_mg_levels 3 -mg_coarse_pc_factor_shift_type nonzero -ksp_view ! ! test: ! suffix: 2 ! nsize: 2 ! args: -ksp_monitor_short -da_grid_x 50 -da_grid_y 50 -pc_type ksp -ksp_ksp_type cg -ksp_pc_type bjacobi -ksp_ksp_rtol 1e-1 -ksp_ksp_monitor -ksp_type pipefgmres -ksp_gmres_restart 5 ! ! test: ! suffix: hyprestruct ! nsize: 3 ! requires: hypre ! args: -ksp_type gmres -pc_type pfmg -dm_mat_type hyprestruct -ksp_monitor -da_refine 3 ! !TEST*/