# 1 "dr_gmres_test.F" ************************************************************************* ** TEST PROGRAM FOR various types of sparse solvers ************************************************************************* ! ! Solving A x = b ! Will also denote x with Phi. ! # 1 "complex_real.h" 1 ! complex_real.h: compile real or complex versions from the same source: # 10 "dr_gmres_test.F" 2 module gm_data ! real*8 :: V_morgan = 1.0 ! real*8 :: V_morgan = 41.0 real*8 :: V_morgan = 0.0 ! nx = # of grid points in 1 direction, including 2 bdy points. ! integer, parameter :: nx = 42 ! integer, parameter :: nx = 59 ! made problems for Neumann once integer, parameter :: nx = 40 integer, parameter :: n = nx**2 real*8 :: dx ! grid spacing real*8 :: xgrid(nx) ! spatial grid complex*16 Phi(nx, nx) complex*16 Phi1d(n) ! equivalent 1-D version of Phi equivalence (Phi,Phi1d) complex*16 Phi_true(nx, nx) ! True solution complex*16 Phi_true1d(n) ! equivalent 1-D version equivalence (Phi_true, Phi_true1d) complex*16 b(nx, nx) ! rhs complex*16 b1d(n) ! equivalent 1-D version of b equivalence (b,b1d) logical :: bc_grad = .true. ! = .false. b.c. on Phi (Dirichlet b.c.'s) ! = .true. b.c. on grad(Phi) (Neumann b.c.'s) real*8 :: pi end module program validation * use gm_data implicit none * integer i, j, m integer info(3) integer nout * complex*16 xtrue(n) real*8 rinfo(2) real*8 x real*8 tol ! relative error tolerance * complex*16 ZERO, ONE parameter (ZERO = (0.0d0, 0.0d0), ONE = (1.0d0, 0.0d0)) ! external Ax_sub_1d ! subroutine to evaluate A.x for 1d Poisson Eq. external Ax_sub_2d ! subroutine to evaluate A.x * logical :: problem_delta = .false. pi = abs(acos(-1.0)) * *************************************************************** ** Generate the test matrix A and set the right-hand side ** in positions (n+1) to 2n of the array work. ** The right-hand side is chosen such that the exact solution ** is the vector of all ones. *************************************************************** * write(*,*) '***********************************************' write(*,*) "This code is an example of an algorithm I'll call" write(*,*) 'optimized dynamic relaxation. Another possible' write(*,*) 'name is recycled GMRES(2), which is like standard' write(*,*) 'GMRES(2) but with a modified restarting method that' write(*,*) 'recycles previous search directions rather than just' write(*,*) & 'throwing them away. This is equivalent to LGMRES(1,1)' write(*,*) '(Loose GMRES), and can also be though of as an' write(*,*) 'optimized 2cd Order Richardson iteration, or' write(*,*) 'a modified form of 2cd order Chebyshev iteration.' write(*,*) 'Results are written in output files' write(*,*) 'fort.* and sol_zTestgmres.' write(*,*) '***********************************************' write(*,*) if(bc_grad) then ! Neumann b.c.s: end grid points lie on either side of ! boundary to give grad(Phi)=0 at bdy. dx = 1.0/(nx-2) do i=1,nx xgrid(i) = -dx/2 + (i-1)*dx enddo else ! Dirichlet b.c.s: end grid points lie on the boundaries: dx = 1.0/(nx-1) do i=1,nx xgrid(i) = (i-1)*dx enddo endif write(*,*) "Matrix size nx=", nx, "N = nx^2 = ", n !---------------------------------------- if(problem_delta) then print *, "can only do Morgan problem at present" stop ! ! or give an explicit b and ignore xtrue: b1d = ZERO ! one or two delta functions: ! b1d(n/2+1) = -1/dx**2 b1d(n*0.75+1) = -1/dx**2 else !---------------------------------------- ! RHS for the Morgan problem: print *, "doing the Morgan problem with V_morgan = ", V_morgan call set_x0_rhs_2d endif !! !----------------------------------- ***************************************** ** 1st or 2cd order Richardson iterative solution of A x = b: ***************************************** * tol = 1.e-9 ! tol = 1.e-14 call dr_gmres(n, Ax_sub_2d, Phi, b, tol, info, rinfo) * ******************************* * dump the solution on a file ******************************* * nout = 11 open(nout,FILE='sol_zTestgmres',STATUS='unknown') ! write(nout,*) 'Restart : ', m write(nout,*) 'info(1) = ',info(1),' info(2) = ',info(2) write(nout,*) 'rinfo(1) = ',rinfo(1),' rinfo(2) = ',rinfo(2) write(nout,*) 'Solution : ' do j=1,n write(nout,*) Phi1d(j) enddo write(nout,*) ' ' * call dump_data write(*,*) "# iterations =", info(2), "relative residual err=", & rinfo(2) write(*,*) "average convergence rate /iteration", & rinfo(2)**(1./info(2)) 100 continue * stop end !----------------------------------------------------------------------- subroutine Ax_sub_1d(n, x, Ax, nx) ! evaluate Matrix-Vector product A x and return the result in Ax implicit none integer n, nx complex*16 x(n), Ax(n) real*8 Delta integer i ! for a 1-D Poisson equation, on the domain [0,1] with a step size 1/(n-1), ! with x=0 boundary conditions. Delta = 1.0/(nx-1) ! skip boundary points at i=1, n do i=2,n-1 Ax(i) = (x(i+1)-2*x(i)+x(i-1))/Delta**2 enddo return end !----------------------------------------------------------------------- subroutine Ax_sub_2d(n, x, Ax) ! evaluate Matrix-Vector product A x and return the result in Ax ! for the 2-D Poisson equation or the isotropic Diffusion problem. use gm_data, only: v_morgan, bc_grad, nx, pi implicit none integer n complex*16 x(nx,nx), Ax(nx,nx) real*8 Delta, unit, dt, D integer i, j ! For 2-D Poisson problem: dt = 1 unit = 0 D = 1 ! ! For isotropic diffusion problem: unit = 1 D = 0.01 ! ! choose time step to scale with dt*D*(pi/L)**2 = 0.2 ! ! i.e., moderately resolve the diffusion rate at the lowest k mode dt = 0.2/(D*pi**2) ! for a 2-D Poisson equation, on the domain [0,1] x [0,1] ! with a step size 1/(n-1), with x=0 boundary conditions. Delta = 1.0/(nx-1) ! skip boundary points at i=1, n do i = 2 , nx-1 do j = 2, nx-1 Ax(i,j) = unit*x(i,j) & - dt*D*((x(i+1,j )-2*x(i,j)+x(i-1,j ))/Delta**2 & +(x(i ,j+1)-2*x(i,j)+x(i ,j-1))/Delta**2 ) & + dt*V_morgan*(x(i+1,j)-x(i-1,j))/(2*Delta) enddo enddo ! set results at boundaries if(.not. bc_grad) then !Dirichlet B.C.'s, (Phi=b) Ax(:,1) = x(:,1) Ax(:,nx) = x(:,nx) Ax(1,:) = x(1,:) Ax(nx,:) = x(nx,:) else ! for a Poisson problem (or dt->infinity limit of diffusion ! problem) with Neumann b.c.'s, explore improving convergence ! by projecting out the null space: ! call clean_null(Ax, Ax_avg) ! call clean_null(Ax, Ax_avg) ! Neumann B.C.'s (grad(Phi)=0 on bdy's) ! Ax(1,2) and Ax(nx,2) are undefined at this point, so the following ! would generate warning errors: ! ! Ax(:,1) = Ax(:,2) ! ! Though this is irrelevant, be careful to use only defined parts of ! the array so that run-time checks for undefined variables ! are useful. Ax(2:nx-1, 1) = Ax(2:nx-1, 2) Ax(2:nx-1, nx) = Ax(2:nx-1, nx-1) Ax(1,:) = Ax(2,:) Ax(nx,:) = Ax(nx-1,:) endif return end subroutine Ax_sub_2d !----------------------------------------------------------------------- subroutine set_x0_rhs_2d ! set initial x and rhs for the 2D Morgan convection-diffusion problem: use gm_data implicit none integer i, j, kmin, kmax, kmid complex*16 phi_avg real*8 x, z1, z2 logical :: random_solution = .true. ! to have a sequence of "random" numbers that is the same ! each time you run the code (which is important to have reproducible ! results that can be debugged), leave the next line commented out. ! To reinitialize the random_seed from the system clock etc., ! uncomment the next line : ! ! call random_seed !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Manufature a random solution, and then find the rhs that ! must correspond to that: do i = 1, nx do j = 1, nx call random_number(z1) call random_number(z2) Phi_true(i,j) = real(cmplx(z1,z2)) ! (above logic works for real or complex versions of Phi_true) enddo enddo ! smooth this solution a few times (makes problem slightly harder): call smooth(Phi_true) call smooth(Phi_true) call smooth(Phi_true) if(bc_grad) then ! well-posed Poisson problem (or dt->infinity limit of diffusion problem) ! with Neumann b.c.'s requires = 0. Set =0 also: call clean_null(Phi_true, Phi_avg) ! impose grad(Phi)=0 b.c.'s: Phi_true(:,1) = Phi_true(:,2) Phi_true(:,nx) = Phi_true(:,nx-1) Phi_true(1,:) = Phi_true(2,:) Phi_true(nx,:) = Phi_true(nx-1,:) endif ! Find the RHS that would correspond to this manufactured solution: call Ax_sub_2d(n, Phi_true1d, b) ! Initial guess: Phi = 0.0 ! Make sure initial guess satisfies b.c.'s: if(.not. bc_grad) then Phi(:,1) = b(:,1) Phi(:,nx) = b(:,nx) Phi(1,:) = b(1,:) Phi(nx,:) = b(nx,:) endif if(random_solution) return !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Other initialization methods: ! ! ! ! initialize the boundary values of b first: ! if(bc_grad) then ! b = 0 ! else ! ! b = 50.0 ! boundary values for x ! b = 0.0 ! boundary values for x ! endif do i = 2, nx-1 do j = 2, nx-1 ! Morgan problem: b(i,j) = -41**2 ! Put in spatially varying temperature for diffusion/Poisson problems: ! with components in the lowest k, highest k, and intermediate k eigenmodes: x = xgrid(i) if(bc_grad) then ! grad(phi)=0 @ bdy's, eigemodes are cos(n*pi*x) for n=0,1,...(nx-1) kmin = 1 kmax = nx-1 kmid = sqrt(1.0*nx) ! If only 2 eigenmodes are excited, then this algorithm converges in ! exactly 2 steps. But if 3 disparate eigenmodes are excited, then ! it converges at a more typical rate. ! 3 exact eigenmodes: Phi_true(i,j) = & cos(kmin*pi*x) / (kmin*pi)**2 & + cos(kmax*pi*x) / (kmax*pi)**2 & - cos(kmid*pi*x) / (kmid*pi)**2 else ! phi=0 @ bdy's, eigemodes are sin(n*pi*x) for n=1,...(nx-2) b(i,j) = sin(pi*x) + sin((nx-2)*pi*x) & + sin(sqrt(1.0*nx)*pi*x) endif enddo enddo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(bc_grad) then ! well-posed Poisson problem (or dt->infinity limit of diffusion problem) ! with Neumann b.c.'s requires = 0. call clean_null(Phi_true, Phi_avg) ! enforce grad(phi)=0 b.c.'s: Phi_true(:,1) = Phi_true(:,2) Phi_true(:,nx) = Phi_true(:,nx-1) Phi_true(1,:) = Phi_true(2,:) Phi_true(nx,:) = Phi_true(nx-1,:) ! Find the RHS that would correspond to this manufactured solution: call Ax_sub_2d(n, Phi_true1d, b) endif ! if(bc_grad) then ! call clean_null(b, b_avg) ! print *, "b_avg =", b_avg ! call clean_null(b, b_avg) ! print *, "b_avg =", b_avg ! endif ! boundary conditions: ! initial x0 must satisfy the b.c.'s too: if(.not. bc_grad) then ! Dirichlet B.C.'s: ! initial guess at solution x0 (or Phi0): Phi = 0 ! Phi = b ! same as Neumann b.c.s for comparison Phi(:,1) = b(:,1) Phi(:,nx) = b(:,nx) Phi(1,:) = b(1,:) Phi(nx,:) = b(nx,:) else ! Neumann B.C.'s (grad(Phi)=0 on bdy's) ! b.c.'s already satisfied by call Ax_sub(n,Phi_true,b) ! b(:,1) = b(:,2) ! b(:,nx) = b(:,nx-1) ! b(1,:) = b(2,:) ! b(nx,:) = b(nx-1,:) Phi = 0 ! Phi = b ! ! Phi(:,1) = Phi(:,2) ! Phi(:,nx) = Phi(:,nx-1) ! Phi(1,:) = Phi(2,:) ! Phi(nx,:) = Phi(nx-1,:) endif return end !----------------------------------------------------------------------- subroutine clean_null(phi, phi_avg) use gm_data, only: nx, bc_grad implicit none complex*16 phi(nx,nx) complex*16 phi_avg integer i, j, i1, j1 phi_avg = 0.0 if(bc_grad) then i1 = 2 j1 = 2 else i1 = 1 j1 = 1 endif do i = i1, nx-1 do j = j1, nx-1 phi_avg = phi_avg + phi(i,j) enddo enddo phi_avg = phi_avg / (nx-i1) / (nx-j1) phi = phi - phi_avg return end subroutine clean_null !----------------------------------------------------------------------- subroutine smooth(phi) use gm_data, only: nx, bc_grad implicit none complex*16 phi(nx,nx) integer i, j, i1, j1 do i = 2, nx-1 do j = 2, nx-1 phi(i,j) = (phi(i,j) & + phi(i+1,j) + phi(i-1,j) & + phi(i,j+1) + phi(i,j-1))/5.0 enddo enddo return end subroutine smooth !----------------------------------------------------------------------- subroutine riter2(n, Ax_sub, x, b, info, rinfo, nx) ! solve A x = b using 2cd order Richardson iteration ! on input, x contains the initial guess. implicit none integer n, nx complex*16 x(n), b(n) integer info(2) real*8 rinfo(2) complex*16 Ax(n), xp(n), xn(n), res(n) complex*16, parameter :: ONE = (1.0d0, 0.0d0) complex*16, parameter :: ZERO = (0.0d0, 0.0d0) real*8 dt, nu, err, pi, kmin, alpha, beta1, nudt_half integer iter, niter, maxiter external Ax_sub dt=0.5*2.0 nu=2*2*3.14159/n * 1.1 maxiter = 1000 pi = 2*abs(acos(0.0)) ! Calculated optimal iteration parameters: kmin=2*pi/(n+1) alpha = 1-2*kmin beta1 = 1/4.+alpha/4+sqrt(alpha)/2 nudt_half=(1-alpha)/(1+alpha) dt = sqrt(beta1*(1+nudt_half)) ! dt = (1+2*kmin)/(1+kmin)/2*(1-kmin+sqrt(1-2*kmin)) ! dt = sqrt(dt) nu=nudt_half*2/dt*0.33 ! must be an error in the above calculation of the optimal parameters, because ! I empirically find a better fit is obtained if I multiply this nu by 0.33... write(*,*) "alpha, beta1 =", alpha, beta1 write(*,*) "dt, nu =", dt, nu write(*,*) "predicted error reduction / iteration =", 1/(1+nu*dt) write(*,*) "b=", b write(15,*) "# iteration #, err" do iter = 1, maxiter ! perform the matrix vector product Ax = A * x call Ax_sub(n, x, Ax, nx) res = Ax - b ! write(*,*) "x =", x ! write(*,*) "res =", res xn = x + (x-xp)*(1-nu*dt/2)/(1+nu*dt/2) + dt**2*res/(1+nu*dt/2) xp = x x = xn err=maxval(abs(res)) write(15,*) iter, err if(err .le. 1.e-10) goto 20 enddo 20 niter = iter info(2)=niter rinfo(2)=err return end subroutine riter2 subroutine dump_data use gm_data implicit none integer i real*8 err2_final, Phi_true2_mag complex*16 zdot ! call clean_null(Phi,Phi_avg) Phi_true2_mag = zdot(Phi_true1d, Phi_true1d, n) err2_final = zdot(Phi-Phi_true, Phi-Phi_true, n) err2_final = sqrt(abs(err2_final)) print *, "L2 relative error |Phi-Phi_true|/|Phi_true| =", & sqrt(err2_final/Phi_true2_mag) write(31,*) "# x Phi(x,j) at j=1" do i = 1, nx write(31,*) xgrid(i), Phi(i,1) enddo write(32,*) "# x Phi(x,j) at j=2" do i = 1, nx write(32,*) xgrid(i), Phi(i,2) enddo write(33,*) "# x Phi(x,j) at j=nx/2, where nx=", nx do i = 1, nx write(33,*) xgrid(i), Phi(i,nx/2) enddo write(34,*) "# x Phi(x,j) at j=nx-1, where nx=", nx do i = 1, nx write(34,*) xgrid(i), Phi(i,nx-1) enddo write(35,*) "# x Phi(x,j) at j=nx, where nx=", nx do i = 1, nx write(35,*) xgrid(i), Phi(i,nx) enddo write(41,*) "# x Phi(x,j) at j=1" do i = 1, nx write(41,*) xgrid(i), Phi_true(i,1) enddo write(42,*) "# x Phi(x,j) at j=2" do i = 1, nx write(42,*) xgrid(i), Phi_true(i,2) enddo write(43,*) "# x Phi(x,j) at j=nx/2, where nx=", nx do i = 1, nx write(43,*) xgrid(i), Phi_true(i,nx/2) enddo write(44,*) "# x Phi(x,j) at j=nx-1, where nx=", nx do i = 1, nx write(44,*) xgrid(i), Phi_true(i,nx-1) enddo write(45,*) "# x Phi(x,j) at j=nx, where nx=", nx do i = 1, nx write(45,*) xgrid(i), Phi_true(i,nx) enddo write(12,*) "# x Phi(x), Phi_true(x)" write(12,*) "0 0.0 0.0" do i=1,n write(12,*) i, real(Phi1d(i)) enddo write(12,*) n+1," 0.0 0.0" return end subroutine dump_data