MODULE green #include 'include.f90' #include "finclude/petscdef.h" USE types USE grid USE petscsys USE petscda USE petscis USE petscmg USE petscvec USE petscmat IMPLICIT NONE CONTAINS !-------------------------------------------- ! subroutine MatrixAntiplaneFV ! computes the forward model of the Poisson's ! equation using the finite volume equations. ! ! sylvain barbot (04/18/11) - original form !-------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "matrixantiplanefv" SUBROUTINE matrixantiplanefv(A,x,y,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: x Vec, INTENT(OUT) :: y PetscErrorCode, INTENT(OUT) :: ierr Vec :: lx,ly PetscInt :: istart,iend,m,n PetscScalar, POINTER :: u(:) PetscScalar, POINTER :: v(:) INTEGER :: rank,isize,its,sw INTEGER :: i2,i3,i2i,i3i,i, & i000,i0p0,i00p, & i0m0,i00m PetscScalar :: dx2l,dx2m,dx2r,dx3l,dx3m,dx3r PetscScalar :: lb PetscScalar :: g,g0m0,g0p0,g00m,g00p PetscScalar :: rtol PetscReal :: rnorm TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) !PRINT *, rank, isize CALL KSPGetTolerances(c%ksp,rtol,PETSC_NULL,PETSC_NULL,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); CALL KSPGetIterationNumber(c%ksp,its,ierr);CHKERRQ(ierr); CALL KSPGetResidualNorm(c%ksp,rnorm,ierr);CHKERRQ(ierr); IF (rank.EQ.0.AND.its.EQ.0) norm_0=rnorm IF (rank.EQ.0.AND.MOD(its,31).EQ.0) PRINT '("iteration ",I5.5, & ", relative residual norm ",ES9.4E1,"/",ES9.4E1)', its,rnorm/norm_0,rtol ! retrieve matrix dimension CALL MatGetSize(A,m,n,ierr);CHKERRQ(ierr); CALL VecGetOwnershipRange(x,istart,iend,ierr);CHKERRQ(ierr); CALL DACreateLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); CALL VecDuplicate(lx,ly,ierr);CHKERRQ(ierr); CALL DAGlobalToLocalBegin(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL DAGlobalToLocalEnd(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(lx,u,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(ly,v,ierr);CHKERRQ(ierr); CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! retrieve stencil width (ghost-point padding) CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & sw,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); i2i=info%xs i3i=info%ys-1 ! balance between Neumann and Dirichlet boundary conditions DO i=istart,iend-1,info%dof i3=(i-istart)/(info%xm*info%dof) i2=(i-istart-i3*info%xm*info%dof)/info%dof i3=i3+i3i ! i3 in (-1 .. sx3-2 ) i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) i000=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof IF (i2.EQ.0) THEN u(i0m0+IU1)=u(i0p0+IU1)*in%phi END IF IF (i2.EQ.in%sx2-1) THEN u(i0p0+IU1)=u(i0m0+IU1)*in%phi END IF IF (i3.EQ.(in%sx3-2)) THEN u(i00p+IU1)=0.d0 END IF END DO ! homogenization constant for boundary condition equation lb=1.d0/MIN(MINVAL(in%dx2),MINVAL(in%dx3)) ! forward model DO i=istart,iend-1,info%dof i3=(i-istart)/(info%xm*info%dof) i2=(i-istart-i3*info%xm*info%dof)/info%dof i3=i3+i3i ! i3 in (-1 .. sx3-2 ) i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) i000=1+(((1+i3-i3i+0)*(info%xm+2)+1+i2-i2i+0))*info%dof i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof i00m=1+(((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof dx2r=in%dx2(i2+1) IF (i2.GT.0) THEN dx2l=in%dx2(i2+1-1) ELSE ! ghost point dx2l=dx2r END IF dx2m=(dx2l+dx2r)/2.d0 dx3r=in%dx3(MAX(0,i3)+1) IF (i3.GT.0) THEN dx3l=in%dx3(i3+1-1) ELSE ! ghost point dx3l=dx3r END IF dx3m=(dx3l+dx3r)/2.d0 ! local shear modulus & Lamé parameter g=moduli(1+(((1+MAX(0,i3-i3i-1))*(info%xm+2)+1+i2-i2i))*NMOD) g0p0=(moduli(1+(((sw+MAX(0,i3-i3i-1))*(info%xm+2*sw)+sw+i2-i2i+1))*NMOD)+g)/2.d0 g0m0=(moduli(1+(((sw+MAX(0,i3-i3i-1))*(info%xm+2*sw)+sw+i2-i2i-1))*NMOD)+g)/2.d0 g00p=(moduli(1+(((sw+MAX(0,i3-i3i-1)+1)*(info%xm+2*sw)+sw+i2-i2i))*NMOD)+g)/2.d0 g00m=(moduli(1+(((sw+MAX(0,i3-i3i-1)-1)*(info%xm+2*sw)+sw+i2-i2i))*NMOD)+g)/2.d0 ! top boundary condition SELECT CASE (i3) CASE (-1) ! do not use the ghost points with finite-volume method v(i000+IU1)=1.d0 CASE (0) ! free surface boundary condition g u1,3 = 0 ! using a control volume starting at the surface v(i000+IU1)=dx3r*(+g0p0/dx2r*( u(i0p0+IU1)-u(i000+IU1) ) & -g0m0/dx2l*( u(i000+IU1)-u(i0m0+IU1) )) & +dx2m*(+g00p/dx3r*( u(i00p+IU1)-u(i000+IU1) )) !v(i000+IU1)=g*(u(i00p+IU1)-u(i000+IU1)) CASE DEFAULT ! interior points in x1-direction, finite-volume equation: ! ! |r2 |r3 ! A2 ( g u1,2 )| + A3 ( g u1,3 )| ! |l2 |l3 ! v(i000+IU1)=dx3m*(+g0p0/dx2r*( u(i0p0+IU1)-u(i000+IU1) ) & -g0m0/dx2l*( u(i000+IU1)-u(i0m0+IU1) )) & +dx2m*(+g00p/dx3r*( u(i00p+IU1)-u(i000+IU1) ) & -g00m/dx3l*( u(i000+IU1)-u(i00m+IU1) )) END SELECT END DO CALL DALocalToGlobal(c%dau,ly,INSERT_VALUES,y,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(lx,u,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(ly,v,ierr);CHKERRQ(ierr); CALL DARestoreLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); CALL DARestoreLocalVector(c%dau,ly,ierr);CHKERRQ(ierr); !CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); !CALL VecDuplicate(x,y,ierr);CHKERRQ(ierr); END SUBROUTINE matrixantiplanefv !-------------------------------------------- ! subroutine MatrixAntiplaneCoarseFV ! computes the forward model of the Poisson's ! equation using the finite volume equations. ! ! sylvain barbot (04/18/11) - original form !-------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "matrixantiplanecoarsefv" SUBROUTINE matrixantiplanecoarsefv(A,x,y,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: x Vec, INTENT(OUT) :: y PetscErrorCode, INTENT(OUT) :: ierr !Vec :: lx,ly PetscInt :: istart,iend,m,n PetscScalar, POINTER :: u(:) PetscScalar, POINTER :: v(:) INTEGER :: rank,isize,i !INTEGER :: its,sw !INTEGER :: i2,i3,i2i,i3i,i, & ! i000,i0p0,i00p, & ! i0m0,i00m !PetscScalar :: dx2l,dx2m,dx2r,dx3l,dx3m,dx3r !PetscScalar :: g,g0m0,g0p0,g00m,g00p !PetscScalar :: rtol !PetscReal :: rnorm !TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) !PRINT *, rank, isize PRINT *, rank, "matrixantiplanecoarsefv" !CALL KSPGetTolerances(c%cksp,rtol,PETSC_NULL,PETSC_NULL,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !CALL KSPGetIterationNumber(c%cksp,its,ierr);CHKERRQ(ierr); !CALL KSPGetResidualNorm(c%cksp,rnorm,ierr);CHKERRQ(ierr); !IF (rank.EQ.0.AND.its.EQ.0) norm_0=rnorm !IF (rank.EQ.0.AND.MOD(its,31).EQ.0) PRINT '("iteration ",I5.5, & ! ", relative residual norm ",ES9.4E1,"/",ES9.4E1)', its,rnorm/norm_0,rtol ! retrieve matrix dimension CALL MatGetSize(A,m,n,ierr);CHKERRQ(ierr); CALL VecGetOwnershipRange(x,istart,iend,ierr);CHKERRQ(ierr); !CALL DACreateLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL VecDuplicate(lx,ly,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalBegin(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalEnd(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(y,v,ierr);CHKERRQ(ierr); !CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! retrieve stencil width (ghost-point padding) !CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! sw,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !i2i=info%xs !i3i=info%ys-1 ! balance between Neumann and Dirichlet boundary conditions !DO i=istart,iend-1,info%dof ! i3=(i-istart)/(info%xm*info%dof) ! i2=(i-istart-i3*info%xm*info%dof)/info%dof ! i3=i3+i3i ! i3 in (-1 .. sx3-2 ) ! i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) ! i000=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof ! i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof ! i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! IF (i2.EQ.0) THEN ! u(i0m0+IU1)=u(i0p0+IU1)*in%phi ! END IF ! IF (i2.EQ.in%sx2-1) THEN ! u(i0p0+IU1)=u(i0m0+IU1)*in%phi ! END IF ! IF (i3.EQ.(in%sx3-2)) THEN ! u(i00p+IU1)=0.d0 ! END IF !END DO ! forward model !DO i=istart,iend-1,info%dof DO i=istart,iend-1,1 ! i3=(i-istart)/(info%xm*info%dof) ! i2=(i-istart-i3*info%xm*info%dof)/info%dof ! i3=i3+i3i ! i3 in (-1 .. sx3-2 ) ! i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) ! i000=1+(((1+i3-i3i+0)*(info%xm+2)+1+i2-i2i+0))*info%dof ! i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof ! i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof ! i00m=1+(((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! dx2r=in%dx2(i2+1) ! IF (i2.GT.0) THEN ! dx2l=in%dx2(i2+1-1) ! ELSE ! ! ghost point ! dx2l=dx2r ! END IF ! dx2m=(dx2l+dx2r)/2.d0 ! dx3r=in%dx3(MAX(0,i3)+1) ! IF (i3.GT.0) THEN ! dx3l=in%dx3(i3+1-1) ! ELSE ! ! ghost point ! dx3l=dx3r ! END IF ! dx3m=(dx3l+dx3r)/2.d0 ! ! local shear modulus & Lamé parameter ! g=moduli(1+(((1+MAX(0,i3-i3i-1))*(info%xm+2)+1+i2-i2i))*NMOD) ! g0p0=(moduli(1+(((sw+MAX(0,i3-i3i-1))*(info%xm+2*sw)+sw+i2-i2i+1))*NMOD)+g)/2.d0 ! g0m0=(moduli(1+(((sw+MAX(0,i3-i3i-1))*(info%xm+2*sw)+sw+i2-i2i-1))*NMOD)+g)/2.d0 ! g00p=(moduli(1+(((sw+MAX(0,i3-i3i-1)+1)*(info%xm+2*sw)+sw+i2-i2i))*NMOD)+g)/2.d0 ! g00m=(moduli(1+(((sw+MAX(0,i3-i3i-1)-1)*(info%xm+2*sw)+sw+i2-i2i))*NMOD)+g)/2.d0 ! top boundary condition ! SELECT CASE (i3) ! CASE (-1) ! ! do not use the ghost points with finite-volume method ! v(i000+IU1)=1.d0 ! CASE (0) ! ! free surface boundary condition g u1,3 = 0 ! ! using a control volume starting at the surface ! v(i000+IU1)=dx3r*(+g0p0/dx2r*( u(i0p0+IU1)-u(i000+IU1) ) & ! -g0m0/dx2l*( u(i000+IU1)-u(i0m0+IU1) )) & ! +dx2m*(+g00p/dx3r*( u(i00p+IU1)-u(i000+IU1) )) ! !v(i000+IU1)=g*(u(i00p+IU1)-u(i000+IU1)) ! CASE DEFAULT ! interior points in x1-direction, finite-volume equation: ! ! |r2 |r3 ! A2 ( g u1,2 )| + A3 ( g u1,3 )| ! |l2 |l3 ! ! v(i000+IU1)=dx3m*(+g0p0/dx2r*( u(i0p0+IU1)-u(i000+IU1) ) & ! -g0m0/dx2l*( u(i000+IU1)-u(i0m0+IU1) )) & ! +dx2m*(+g00p/dx3r*( u(i00p+IU1)-u(i000+IU1) ) & ! -g00m/dx3l*( u(i000+IU1)-u(i00m+IU1) )) ! END SELECT v(i-istart+1)=u(i-istart+1) END DO !CALL DALocalToGlobal(c%dau,ly,INSERT_VALUES,y,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(y,v,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,ly,ierr);CHKERRQ(ierr); !CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); !CALL VecDuplicate(x,y,ierr);CHKERRQ(ierr); END SUBROUTINE matrixantiplanecoarsefv !--------------------------------------------------------------------- ! subroutine greenfunction ! solves the elasto-static problem A*x = b with previously-defined ! matrx A and vector b. ! ! sylvain barbot (08-20-10) - original form !--------------------------------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "greenfunction" SUBROUTINE greenfunction(c) TYPE(COMP_STRUC), INTENT(INOUT) :: c INTEGER :: ierr CALL KSPSolve(c%ksp,c%b(1),c%x(1),ierr);CHKERRQ(ierr); END SUBROUTINE greenfunction !--------------------------------------------------------------------- !> subroutine green_init_mg !! creates the matrix for the 2-D Navier elasto-static equation for !! multi-grid preconditioner. !! !! per petsc conventions, !! !! 0 is always the coarse grid !! !! n-1 is always the fine grid !! !! there is no restriction from level 0 to -1. !! !! attention: !! !! MGGetGoarseSolve(pc,KSP *cksp) == MGGetSmoother(pc,0,KSP *cksp) !! !! There are 3 work vectors for multi grid. b and X are not needed !! on finest grid. r is not needed on coarsest grid. ! ! sylvain barbot (05-03-11) - original form !--------------------------------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "green_init_mg" SUBROUTINE green_init_mg(c) TYPE(COMP_STRUC), INTENT(INOUT) :: c PetscMPIInt :: rank,isize,ierr PetscReal :: rtol,abstol,dtol PetscInt :: dummy,x,y,z,maxits PetscInt :: pm,pn INTEGER :: l,p,iostatus,count PetscInt, DIMENSION(:), ALLOCATABLE :: lx,ly TYPE(DALOCALINFOF90) :: info,cinfo,minfo KSP :: sksp CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) ! setting the number of multi-grid levels c%nlevels=2 !------------------------------------- ! higher-level multi-grid solver !------------------------------------- CALL KSPCreate(PETSC_COMM_WORLD,c%ksp,ierr);CHKERRQ(ierr); CALL KSPGetPC(c%ksp,c%pc,ierr);CHKERRQ(ierr); CALL PCSetType(c%pc,PCMG,ierr);CHKERRQ(ierr); CALL PCMGSetLevels(c%pc,c%nlevels,PETSC_NULL_OBJECT,ierr);CHKERRQ(ierr); ! change with option -pc_mg_type [multiplicative, additive, full, kaskade] CALL PCMGSetType(c%pc,PC_MG_MULTIPLICATIVE,ierr);CHKERRQ(ierr); ! change with option -pc_mg_cycle_type v or w ! cycle types are PC_MG_CYCLE_V or PC_MG_CYCLE_W CALL PCMGSetCycleType(c%pc,PC_MG_CYCLE_V,ierr);CHKERRQ(ierr); ! control the amount of pre- and postsmoothing with the options ! -pc_mg_smoothup m and -pc_mg_smoothdown n CALL PCMGSetNumberSmoothUp(c%pc,1,ierr);CHKERRQ(ierr); CALL PCMGSetNumberSmoothDown(c%pc,1,ierr);CHKERRQ(ierr); ! dimension of finest grid CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! find out number of processors per direction CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & pm,pn,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL,PETSC_NULL,ierr);CHKERRQ(ierr); !-------------------------------------------------------------------- ! coarse-level direct solver !-------------------------------------------------------------------- CALL KSPCreate(PETSC_COMM_WORLD,c%cksp,ierr);CHKERRQ(ierr); ! dummy is info%xm CALL DAGetCorners(c%dau,x,y,PETSC_NULL_INTEGER, & dummy,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); ! global size at coarsest level cinfo%mx=info%mx/2**(c%nlevels-1)+1 cinfo%my=info%my/2**(c%nlevels-1)+1 cinfo%mz=1 ! locally owned size at coarsest level CALL getmxl(x,info%xm,c%nlevels-1,cinfo%xm) CALL getmxl(y,info%ym,c%nlevels-1,cinfo%ym) cinfo%zm=1 cinfo%dof=info%dof ! square matrix CALL MatCreateShell(PETSC_COMM_WORLD, & cinfo%zm*cinfo%ym*cinfo%xm*info%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof, & cinfo%mz*cinfo%my*cinfo%mx*info%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof, & PETSC_NULL,c%cA,ierr);CHKERRQ(ierr); CALL MatShellSetOperation(c%cA,MATOP_MULT,matrixantiplanecoarsefv,ierr);CHKERRQ(ierr); CALL MatAssemblyBegin(c%cA,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL MatAssemblyEnd(c%cA,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL KSPSetOperators(c%cksp,c%cA,c%cA,DIFFERENT_NONZERO_PATTERN,ierr);CHKERRQ(ierr); CALL KSPSetUp(c%cksp,ierr);CHKERRQ(ierr); CALL PCMGGetCoarseSolve(c%pc,c%cksp,ierr);CHKERRQ(ierr); !-------------------------------------------------------------------- ! interpolation and restriction for all but the coarsest level !-------------------------------------------------------------------- ! allocate interpolation and restriction matrices ALLOCATE(c%interp(c%nlevels-1),c%restrict(c%nlevels-1),STAT=iostatus) IF (iostatus>0) STOP "could not allocate memory to make interpolation and restriction matrices" DO l=1,c%nlevels-1 ! global size at current level cinfo%mx=(info%mx/2**(c%nlevels-l))*2+1 cinfo%my=(info%my/2**(c%nlevels-l))*2+1 cinfo%mz=1 ! locally-owned size at current level CALL getmxl(x,info%xm,c%nlevels-l-1,cinfo%xm) CALL getmxl(y,info%ym,c%nlevels-l-1,cinfo%ym) cinfo%zm=1 cinfo%dof=info%dof ! global size at previous, coarser, level minfo%mx=info%mx/2**(c%nlevels-l)+1 minfo%my=info%my/2**(c%nlevels-l)+1 minfo%mz=1 ! locally-owned size at previous, coarser, level CALL getmxl(x,info%xm,c%nlevels-l,minfo%xm) CALL getmxl(y,info%ym,c%nlevels-l,minfo%ym) minfo%zm=1 minfo%dof=info%dof ! shell matrix MxN with M>N CALL MatCreateShell(PETSC_COMM_WORLD, & cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,minfo%zm*minfo%ym*minfo%xm*minfo%dof, & cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,minfo%mz*minfo%my*minfo%mx*minfo%dof, & PETSC_NULL,c%interp(l),ierr);CHKERRQ(ierr); CALL MatShellSetOperation(c%interp(l),MATOP_MULT,matrixinterp,ierr);CHKERRQ(ierr); CALL MatAssemblyBegin(c%interp(l),MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL MatAssemblyEnd(c%interp(l),MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); ! shell matrix MxN with M0) STOP "could not allocate memory to make smoothing matrices" ! allocate forward model matrices ALLOCATE(c%rA(c%nlevels),STAT=iostatus) IF (iostatus>0) STOP "could not allocate memory to make forward-model matrices" DO l=0,c%nlevels-1 ! global size at current level cinfo%mx=(info%mx/2**(c%nlevels-l))*2+1 cinfo%my=(info%my/2**(c%nlevels-l))*2+1 cinfo%mz=1 ! locally-owned size at current level CALL getmxl(x,info%xm,c%nlevels-l-1,cinfo%xm) CALL getmxl(y,info%ym,c%nlevels-l-1,cinfo%ym) cinfo%zm=1 cinfo%dof=info%dof !----------------------------------------------------------------------------- ! smoothing matrix !----------------------------------------------------------------------------- ! square smoothing matrix CALL MatCreateShell(PETSC_COMM_WORLD, & cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof, & cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof, & PETSC_NULL,c%sA(l+1),ierr);CHKERRQ(ierr); CALL MatShellSetOperation(c%sA(l+1),MATOP_SOR,matrixsmoothfv,ierr);CHKERRQ(ierr); CALL MatAssemblyBegin(c%sA(l+1),MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL MatAssemblyEnd(c%sA(l+1),MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL PCMGGetSmoother(c%pc,l,sksp,ierr);CHKERRQ(ierr); CALL KSPCreate(PETSC_COMM_WORLD,sksp,ierr);CHKERRQ(ierr); CALL KSPSetOperators(sksp,c%sA(l+1),c%sA(l+1),DIFFERENT_NONZERO_PATTERN,ierr);CHKERRQ(ierr); CALL KSPSetUp(sksp,ierr);CHKERRQ(ierr); !----------------------------------------------------------------------------- ! residual operator !----------------------------------------------------------------------------- ! square forward-model matrix CALL MatCreateShell(PETSC_COMM_WORLD, & cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof,cinfo%zm*cinfo%ym*cinfo%xm*cinfo%dof, & cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof,cinfo%mz*cinfo%my*cinfo%mx*cinfo%dof, & PETSC_NULL,c%rA(l+1),ierr);CHKERRQ(ierr); CALL MatShellSetOperation(c%rA(l+1),MATOP_SOR,matrixresidualfv,ierr);CHKERRQ(ierr); CALL MatAssemblyBegin(c%rA(l+1),MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL MatAssemblyEnd(c%rA(l+1),MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); ! compute formula b - Ax with matrix-free matrix ! PCMGDefaultResidual(Mat mat,Vec b,Vec x,Vec r) !CALL PCMGSetResidual(c%pc,l,PCMGDefaultResidual,c%rA(l+1),ierr);CHKERRQ(ierr); CALL PCMGSetResidual(c%pc,l,matrixresidualfv,c%rA(l+1),ierr);CHKERRQ(ierr); END DO !------------------------------------------------------------- ! DAs for all levels !------------------------------------------------------------- ! allocate work space ALLOCATE(c%dal(c%nlevels),lx(pm),ly(pn),STAT=iostatus) IF (iostatus>0) STOP "could not allocate memory to make DAs." DO l=c%nlevels-1,0,-1 ! global size at current level cinfo%mx=(info%mx/2**(c%nlevels-l+1))*2+1 cinfo%my=(info%my/2**(c%nlevels-l+1))*2+1 cinfo%mz=1 CALL getdatalayout(info%xm,x,y,c%nlevels-l,pm,lx) CALL getdatalayout(info%ym,y,x,c%nlevels-l,pn,ly) CALL DACreate2d(PETSC_COMM_WORLD,DA_XYPERIODIC,DA_STENCIL_BOX,cinfo%mx,cinfo%my, & pm,pn,DGF,1, & lx,ly,c%dal(l+1),ierr);CHKERRQ(ierr); END DO DEALLOCATE(lx,ly) !------------------------------------------------------------- ! work space for all levels ! b and X are not needed on finest grid (l=n-1) ! r is not needed on coarsest grid (l=0) !------------------------------------------------------------- ! allocate work space ALLOCATE(c%b(c%nlevels),c%x(c%nlevels),c%r(c%nlevels),STAT=iostatus) IF (iostatus>0) STOP "could not allocate memory to make work space" DO l=c%nlevels-1,0,-1 IF (l.NE.0) THEN ! create r vector CALL DACreateGlobalVector(c%dal(l+1),c%r(l+1),ierr);CHKERRQ(ierr); CALL PCMGSetR(c%pc,l,c%r(l+1),ierr);CHKERRQ(ierr); END IF IF (l.NE.c%nlevels-1) THEN ! create b and x vectors CALL DACreateGlobalVector(c%dal(l+1),c%b(l+1),ierr);CHKERRQ(ierr); CALL DACreateGlobalVector(c%dal(l+1),c%x(l+1),ierr);CHKERRQ(ierr); ! PC references these vectors so call VecDestroy() when you are finished with them CALL PCMGSetRhs(c%pc,l,c%b(l+1),ierr);CHKERRQ(ierr); CALL PCMGSetX(c%pc,l,c%x(l+1),ierr);CHKERRQ(ierr); END IF END DO ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CALL MatCreateShell(PETSC_COMM_WORLD, & info%zm*info%ym*info%xm*info%dof,info%zm*info%ym*info%xm*info%dof, & info%mz*info%my*info%mx*info%dof,info%mz*info%my*info%mx*info%dof, & PETSC_NULL,c%A,ierr);CHKERRQ(ierr); CALL MatShellSetOperation(c%A,MATOP_MULT,matrixantiplanefv,ierr);CHKERRQ(ierr); CALL MatAssemblyBegin(c%A,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL MatAssemblyEnd(c%A,MAT_FINAL_ASSEMBLY,ierr);CHKERRQ(ierr); CALL KSPSetOperators(c%ksp,c%A,c%A,DIFFERENT_NONZERO_PATTERN,ierr);CHKERRQ(ierr); !CALL KSPGetTolerances(c%ksp,rtol,abstol,dtol,maxits,ierr);CHKERRQ(ierr); !rtol=1.d-4 !CALL KSPSetTolerances(c%ksp,rtol,abstol,dtol,maxits,ierr);CHKERRQ(ierr); ! use -ksp_type bcgs, lgmres, gmres, bcgsl, cgs, tfqmr, cr, minres, lcd, gcr, ngmres CALL KSPSetFromOptions(c%ksp,ierr);CHKERRQ(ierr); CALL KSPSetUp(c%ksp,ierr);CHKERRQ(ierr); !CALL KSPView(c%ksp,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); PRINT *, rank, "done with init green" END SUBROUTINE green_init_mg !-------------------------------------------- ! subroutine MatrixInterp ! interpolates solution from a coarse level ! to a finer level. ! ! sylvain barbot (05/02/11) - original form !-------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "matrixinterp" SUBROUTINE matrixinterp(A,x,y,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: x Vec, INTENT(OUT) :: y PetscErrorCode, INTENT(OUT) :: ierr !Vec :: lx,ly PetscInt :: istart,iend,m,n PetscScalar, POINTER :: u(:) PetscScalar, POINTER :: v(:) INTEGER :: rank,isize,i !INTEGER :: its,sw !INTEGER :: i2,i3,i2i,i3i, & ! i000,i0p0,i00p, & ! i0m0,i00m !PetscScalar :: rtol !PetscReal :: rnorm !TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) !PRINT *, rank, isize PRINT *, rank, "matrixinterp" !CALL KSPGetTolerances(c%ksp,rtol,PETSC_NULL,PETSC_NULL,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !CALL KSPGetIterationNumber(c%ksp,its,ierr);CHKERRQ(ierr); !CALL KSPGetResidualNorm(c%ksp,rnorm,ierr);CHKERRQ(ierr); !IF (rank.EQ.0.AND.its.EQ.0) norm_0=rnorm !IF (rank.EQ.0.AND.MOD(its,31).EQ.0) PRINT '("iteration ",I5.5, & ! ", relative residual norm ",ES9.4E1,"/",ES9.4E1)', its,rnorm/norm_0,rtol ! retrieve matrix dimension CALL MatGetSize(A,m,n,ierr);CHKERRQ(ierr); CALL VecGetOwnershipRange(y,istart,iend,ierr);CHKERRQ(ierr); !CALL DACreateLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL VecDuplicate(lx,ly,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalBegin(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalEnd(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(y,v,ierr);CHKERRQ(ierr); !CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! retrieve stencil width (ghost-point padding) !CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! sw,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !i2i=info%xs !i3i=info%ys-1 ! forward model !DO i=istart,iend-1,info%dof DO i=istart,iend-1,1 ! i3=(i-istart)/(info%xm*info%dof) ! i2=(i-istart-i3*info%xm*info%dof)/info%dof ! i3=i3+i3i ! i3 in (-1 .. sx3-2 ) ! i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) ! i000=1+(((1+i3-i3i+0)*(info%xm+2)+1+i2-i2i+0))*info%dof ! i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof ! i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof ! i00m=1+(((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof v(i-istart+1)=0 END DO !CALL DALocalToGlobal(c%dau,ly,INSERT_VALUES,y,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(y,v,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,ly,ierr);CHKERRQ(ierr); !CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); !CALL VecDuplicate(x,y,ierr);CHKERRQ(ierr); END SUBROUTINE matrixinterp !-------------------------------------------- ! subroutine MatrixRestrict ! restrict (subsample) residuals from a coarse level ! to a finer level. ! ! sylvain barbot (05/02/11) - original form !-------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "matrixrestrict" SUBROUTINE matrixrestrict(A,x,y,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: x Vec, INTENT(OUT) :: y PetscErrorCode, INTENT(OUT) :: ierr !Vec :: lx,ly PetscInt :: istart,iend,m,n PetscScalar, POINTER :: u(:) PetscScalar, POINTER :: v(:) INTEGER :: rank,isize,i !INTEGER :: its,sw !INTEGER :: i2,i3,i2i,i3i, & ! i000,i0p0,i00p, & ! i0m0,i00m !PetscScalar :: rtol !PetscReal :: rnorm !TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) !PRINT *, rank, isize PRINT *, rank, "matrixrestrict" ! retrieve matrix dimension CALL MatGetSize(A,m,n,ierr);CHKERRQ(ierr); CALL VecGetOwnershipRange(y,istart,iend,ierr);CHKERRQ(ierr); !CALL DACreateLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL VecDuplicate(lx,ly,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalBegin(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalEnd(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(y,v,ierr);CHKERRQ(ierr); !CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! retrieve stencil width (ghost-point padding) !CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! sw,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !i2i=info%xs !i3i=info%ys-1 ! forward model !DO i=istart,iend-1,info%dof DO i=istart,iend-1,1 ! i3=(i-istart)/(info%xm*info%dof) ! i2=(i-istart-i3*info%xm*info%dof)/info%dof ! i3=i3+i3i ! i3 in (-1 .. sx3-2 ) ! i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) ! i000=1+(((1+i3-i3i+0)*(info%xm+2)+1+i2-i2i+0))*info%dof ! i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof ! i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof ! i00m=1+(((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof v(i-istart+1)=0 END DO !CALL DALocalToGlobal(c%dau,ly,INSERT_VALUES,y,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(y,v,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,ly,ierr);CHKERRQ(ierr); !CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); !CALL VecDuplicate(x,y,ierr);CHKERRQ(ierr); END SUBROUTINE matrixrestrict !-------------------------------------------------- ! subroutine MatrixSmoothFV ! performed a weighted jacobi scheme based on the ! finite-volume equations of elasticity. ! ! sylvain barbot (05/02/11) - original form !-------------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "matrixsmoothfv" SUBROUTINE matrixsmoothfv(A,b,omega,flag,shift,its,lits,x,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: x Vec, INTENT(OUT) :: b PetscReal, INTENT(IN) :: omega,shift MatSORType, INTENT(IN) :: flag PetscInt, INTENT(IN) :: its, lits PetscErrorCode, INTENT(OUT) :: ierr !Vec :: lx,lb PetscInt :: istart,iend,m,n PetscScalar, POINTER :: u(:) PetscScalar, POINTER :: v(:) INTEGER :: rank,isize,i !INTEGER :: its,sw !INTEGER :: i2,i3,i2i,i3i, & ! i000,i0p0,i00p, & ! i0m0,i00m !PetscScalar :: rtol !PetscReal :: rnorm !TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) !PRINT *, rank, isize PRINT *, rank, "matrixsmoothfv" !CALL KSPGetTolerances(c%ksp,rtol,PETSC_NULL,PETSC_NULL,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !CALL KSPGetIterationNumber(c%ksp,its,ierr);CHKERRQ(ierr); !CALL KSPGetResidualNorm(c%ksp,rnorm,ierr);CHKERRQ(ierr); !IF (rank.EQ.0.AND.its.EQ.0) norm_0=rnorm !IF (rank.EQ.0.AND.MOD(its,31).EQ.0) PRINT '("iteration ",I5.5, & ! ", relative residual norm ",ES9.4E1,"/",ES9.4E1)', its,rnorm/norm_0,rtol ! retrieve matrix dimension CALL MatGetSize(A,m,n,ierr);CHKERRQ(ierr); CALL VecGetOwnershipRange(x,istart,iend,ierr);CHKERRQ(ierr); !CALL DACreateLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL VecDuplicate(lx,lb,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalBegin(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalEnd(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(b,v,ierr);CHKERRQ(ierr); !CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! retrieve stencil width (ghost-point padding) !CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! sw,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !i2i=info%xs !i3i=info%ys-1 ! forward model !DO i=istart,iend-1,info%dof DO i=istart,iend-1,1 ! i3=(i-istart)/(info%xm*info%dof) ! i2=(i-istart-i3*info%xm*info%dof)/info%dof ! i3=i3+i3i ! i3 in (-1 .. sx3-2 ) ! i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) ! i000=1+(((1+i3-i3i+0)*(info%xm+2)+1+i2-i2i+0))*info%dof ! i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof ! i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof ! i00m=1+(((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof v(i-istart+1)=u(i-istart+1) END DO !CALL DALocalToGlobal(c%dau,lb,INSERT_VALUES,b,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(x,u,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(b,v,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,lb,ierr);CHKERRQ(ierr); !CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); !CALL VecDuplicate(x,b,ierr);CHKERRQ(ierr); END SUBROUTINE matrixsmoothfv !-------------------------------------------------- ! subroutine MatrixResidualFV ! evaluates the residual ! ! r = b - A x ! ! using a shell matrix. ! ! complies with the interface ! ! PCMGDefaultResidual(Mat mat,Vec b,Vec x,Vec r) ! ! sylvain barbot (05/02/11) - original form !-------------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "matrixresidualfv" SUBROUTINE matrixresidualfv(A,b,x,r,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: b,x Vec, INTENT(OUT) :: r PetscErrorCode, INTENT(OUT) :: ierr !Vec :: lx,ly PetscInt :: istart,iend,m,n PetscScalar, POINTER :: pb(:) PetscScalar, POINTER :: px(:) PetscScalar, POINTER :: pr(:) INTEGER :: rank,isize,i !INTEGER :: its,sw !INTEGER :: i2,i3,i2i,i3i, & ! i000,i0p0,i00p, & ! i0m0,i00m !PetscScalar :: rtol !PetscReal :: rnorm !TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) !PRINT *, rank, isize PRINT *, rank, "matrixsmoothfv" !CALL KSPGetTolerances(c%ksp,rtol,PETSC_NULL,PETSC_NULL,PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !CALL KSPGetIterationNumber(c%ksp,its,ierr);CHKERRQ(ierr); !CALL KSPGetResidualNorm(c%ksp,rnorm,ierr);CHKERRQ(ierr); !IF (rank.EQ.0.AND.its.EQ.0) norm_0=rnorm !IF (rank.EQ.0.AND.MOD(its,31).EQ.0) PRINT '("iteration ",I5.5, & ! ", relative residual norm ",ES9.4E1,"/",ES9.4E1)', its,rnorm/norm_0,rtol ! retrieve matrix dimension CALL MatGetSize(A,m,n,ierr);CHKERRQ(ierr); CALL VecGetOwnershipRange(b,istart,iend,ierr);CHKERRQ(ierr); !CALL DACreateLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL VecDuplicate(lx,ly,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalBegin(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); !CALL DAGlobalToLocalEnd(c%dau,x,INSERT_VALUES,lx,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(b,pb,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(x,px,ierr);CHKERRQ(ierr); CALL VecGetArrayF90(r,pr,ierr);CHKERRQ(ierr); !CALL DAGetLocalInfoF90(c%dau,info,ierr);CHKERRQ(ierr); ! retrieve stencil width (ghost-point padding) !CALL DAGetInfo(c%dau,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & ! sw,PETSC_NULL_INTEGER, & ! PETSC_NULL_INTEGER,ierr);CHKERRQ(ierr); !i2i=info%xs !i3i=info%ys-1 ! forward model !DO i=istart,iend-1,info%dof DO i=istart,iend-1,1 ! i3=(i-istart)/(info%xm*info%dof) ! i2=(i-istart-i3*info%xm*info%dof)/info%dof ! i3=i3+i3i ! i3 in (-1 .. sx3-2 ) ! i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) ! i000=1+(((1+i3-i3i+0)*(info%xm+2)+1+i2-i2i+0))*info%dof ! i0p0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1))*info%dof ! i00p=1+(((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof ! i0m0=1+(((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1))*info%dof ! i00m=1+(((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0))*info%dof pr(i-istart+1)=pb(i-istart+1)-px(i-istart+1) END DO !CALL DALocalToGlobal(c%dau,ly,INSERT_VALUES,y,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(b,pb,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(x,px,ierr);CHKERRQ(ierr); CALL VecRestoreArrayF90(r,pr,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,lx,ierr);CHKERRQ(ierr); !CALL DARestoreLocalVector(c%dau,ly,ierr);CHKERRQ(ierr); !CALL VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr);CHKERRQ(ierr); !CALL VecDuplicate(x,y,ierr);CHKERRQ(ierr); END SUBROUTINE matrixresidualfv !------------------------------------------------------------------- !> subroutine getDataLayout !! establish the number of elements locally owned by each !! cell at a given level of the multi grid. !! !! allgather corner location of all threads in y direction !! allgather number of locally owned points in x direction !! keep only numbers with unique corner position in y. !! !! @param mx is the number of locally owned elements !! @param cy is the corner position in direction y !! @param l is the multi-grid level !! @param proc is the number of processors in the x direction !! and the size of lx !! @param lx is the output data layout of all threads, !! an integer array with a preallocated size corresponding !! to the number of threads. ! ! sylvain barbot (05/01/11) - original form !------------------------------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "getdatalayout" SUBROUTINE getdatalayout(mx,cx,cy,l,proc,lx) INTEGER, INTENT(IN) :: mx,cx,cy,l,proc INTEGER, DIMENSION(proc), INTENT(OUT) :: lx INTEGER :: rank,isize,iostatus,ierr INTEGER :: i,j ! reference corner position INTEGER :: refcy ! list of all corner positions in the y direction INTEGER, DIMENSION(:), ALLOCATABLE :: allmx,allcx,allcy ! list of relevant corner positions in the x direction INTEGER, DIMENSION(:), ALLOCATABLE :: rcx CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) ALLOCATE(allmx(isize),allcx(isize),allcy(isize),rcx(proc),STAT=iostatus) IF (iostatus>0) STOP "could not allocate corner positions" CALL MPI_ALLGATHER(mx,1,MPI_INTEGER,allmx,1,MPI_INTEGER,PETSC_COMM_WORLD,ierr) CALL MPI_ALLGATHER(cx,1,MPI_INTEGER,allcx,1,MPI_INTEGER,PETSC_COMM_WORLD,ierr) CALL MPI_ALLGATHER(cy,1,MPI_INTEGER,allcy,1,MPI_INTEGER,PETSC_COMM_WORLD,ierr) refcy=allcy(1) j=0 DO i=1,isize IF (allcy(i).EQ.refcy) THEN j=j+1 lx(j)=allmx(i) rcx(j)=allcx(i) END IF END DO IF (proc.NE.j) THEN PRINT *, rank,proc, j STOP "error looping through thread info" END IF ! get the subsampled size at level l DO i=1,proc !PRINT *, "corner=",rcx(i) CALL getmxl(rcx(i),allmx(i),l,lx(i)) END DO DEALLOCATE(allmx,allcx,allcy,rcx) END SUBROUTINE getdatalayout !------------------------------------------------------------------- !> subroutine GetMXL !! evaluates the number of nodes owned by a thread at level !! l of the multi-grid. !! !! @param x is the starting index of locally owned parameters !! @param mx is the size of locally owned data !! @param l is the number of times the array was down-sampled !! @param mxl is the number of locally owned data after l reductions ! ! sylvain barbot (05-04-11) - original form !------------------------------------------------------------------- #undef __FUNCT__ #define __FUNCT__ "getmxl" SUBROUTINE getmxl(x,mx,l,mxl) INTEGER, INTENT(IN) :: x,mx,l INTEGER, INTENT(OUT) :: mxl INTEGER :: k,cx,cmx cx=x cmx=mx DO k=1,l cmx=(cmx-mod(cx,2)+1)/2 cx=(cx+1)/2 END DO mxl=cmx END SUBROUTINE getmxl END MODULE green