program test_relax c--------------------------------------------------------- c Thu Feb 15 10:34:23 PST 2001 c c Test program for routine relax. c For the Black String Coalition. c--------------------------------------------------------- implicit none character*10 cdnm parameter ( cdnm = 'test_relax' ) logical ltrace parameter ( ltrace = .true. ) integer maxdim parameter ( maxdim = 1000 ) real*8 alpha(maxdim, maxdim) real*8 a(maxdim, maxdim) real*8 c(maxdim, maxdim) real*8 br(maxdim, maxdim) real*8 bz(maxdim, maxdim) real*8 krr(maxdim, maxdim) real*8 phi(maxdim, maxdim) real*8 pi(maxdim, maxdim) real*8 phi_r(maxdim, maxdim) real*8 phi_z(maxdim, maxdim) real*8 r(maxdim), z(maxdim) real*8 mg_rhs1(maxdim, maxdim) real*8 mg_rhs2(maxdim, maxdim) real*8 mg_rhs3(maxdim, maxdim) real*8 mg_rhs4(maxdim, maxdim) integer level_r, level_z integer Nr, Nz integer order, maxiter real*8 tol integer i, j integer iargc, i4arg real*8 r8arg real*8 amp, rac, rw real*8 rmax, zmax real*8 dr, dz c========================================================= c========================================================= if ( iargc() .gt. 8 ) goto 999 order = i4arg(1,1) level_r = i4arg(2,5) level_z = i4arg(3,5) amp = r8arg(4,0.0d-1) rac = r8arg(5,5.0d0) rw = r8arg(6,1.0d0) rmax = r8arg(7,1.0d0) zmax = r8arg(8,1.0d0) if ( order .lt. 0 ) goto 999 if ( level_r .lt. 1 ) goto 999 if ( level_z .lt. 1 ) goto 999 Nr = 2**level_r + 1 Nz = 2**level_z + 1 dr = rmax / (Nr-1) dz = zmax / (Nz-1) if ( ltrace ) then write(0,*) cdnm,': order=',order write(0,*) cdnm,': level_r=',level_r write(0,*) cdnm,': level_z=',level_z write(0,*) cdnm,': N_r=',Nr write(0,*) cdnm,': N_z=',Nz write(0,*) cdnm,': amp=',amp write(0,*) cdnm,': rc=',rac write(0,*) cdnm,': rw=',rw write(0,*) cdnm,': rmax=',rmax write(0,*) cdnm,': zmax=',zmax endif call init_geo(alpha,a,c,krr,br,bz,phi,pi, & phi_r,phi_z,Nr,Nz,r,z,rac, & rw,amp,dr, dz, mg_rhs1, & mg_rhs2, mg_rhs3, mg_rhs4) tol = 1.0d-15 maxiter =20000 if ( ltrace ) then write(0,*) cdnm,': maxiter=',maxiter endif call relax(alpha, a, c, br, bz, krr, & phi, pi, phi_r, phi_z, & Nr, Nz, r, z, & order, tol, maxiter, mg_rhs1, & mg_rhs2, mg_rhs3, mg_rhs4) stop 999 continue write(0,*) 'usage: ',cdnm,' '// & ' ' end c##################################################################### c##################################################################### subroutine init_geo(alpha,a,c,krr,br,bz,phi,pi, & phi_r,phi_z,Nr,Nz,r,z,rac, & rw,amp,dr,dz, mg_rhs1, mg_rhs2, & mg_rhs3, mg_rhs4) implicit none integer Nr, Nz real*8 alpha(Nr, Nz) real*8 a(Nr, Nz) real*8 c(Nr, Nz) real*8 krr(Nr, Nz) real*8 br(Nr, Nz) real*8 bz(Nr, Nz) real*8 phi(Nr, Nz) real*8 phi_r(Nr, Nz) real*8 phi_z(Nr, Nz) real*8 pi(Nr, Nz) real*8 r(Nr), z(Nz) real*8 rac, rw, amp integer i, j real*8 dr, dz real*8 cpi parameter ( cpi = 3.141 5926 5358 9793 d0 ) real*8 mg_rhs1(Nr,Nz) real*8 mg_rhs2(Nr,Nz) real*8 mg_rhs3(Nr,Nz) real*8 mg_rhs4(Nr,Nz) c===================================================================== c===================================================================== r(1) = 0.0d0 do i = 2 , Nr r(i) = r(i-1) + dr enddo z(1) = 0.0d0 do j = 2 , Nz z(j) = z(j-1) + dz enddo do i = 1 , Nr do j = 1 , Nz alpha(i, j) = 1.0d0 a(i, j) = 1.0d0 c(i, j) = 1.0d0 krr(i,j) = 0.0d0 br(i,j) = 0.0d0 bz(i,j) = 0.0d0 phi(i,j) = amp * exp(-(r(i)-rac)**2/rw**2) c & * exp(-(z(j)-z(Nz)/2.0d0)**2/rw**2) phi_r(i,j) =-amp * 2.0d0 * (r(i)-rac)/rw**2 * & exp(-(r(i)-rac)**2/rw**2) c & * exp(-(z(j)-z(Nz)/2.0d0)**2/rw**2) phi_z(i,j) = 0.0d0 c phi_z(i,j) = amp * exp(-(r(i)-rac)**2/rw**2) c & * exp(-(z(j)-z(Nz)/2.0d0)**2/rw**2) c & * 2.0d0 * (z(j)-z(Nz)/2.0d0)/rw**2 pi(i,j) = 0.0d0 c pi(i,j) = amp * exp(-(r(i)-rac)**2/rw**2) c & * exp(-(z(j)-z(Nz)/2.0d0)**2/rw**2) mg_rhs1(i,j) =0.0d0 mg_rhs2(i,j) =0.0d0 mg_rhs3(i,j) =0.0d0 mg_rhs4(i,j) =0.0d0 enddo enddo return end