#include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" c note, some function calls in here may not work. pal subroutine linear_ham_solver(CCTK_FARGUMENTS,my_elliptic_tol) implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS INTEGER CCTK_Equals INTEGER, DIMENSION(3) :: sym INTEGER ierr integer icount,i,j,k,bbox_save(6) CCTK_REAL max_resid,one CCTK_REAL my_elliptic_tol integer octant, bitant, quadrant INTEGER reduce_handle, return INTEGER Mlin_index,Nsrc_index,afield_index INTEGER metpsi_index(7) CCTK_REAL AbsTol(3), RelTol(3) INTEGER petsc, sor, jacobi, bamlin, nonlin INTEGER ivp_psi_linear_old_index CCTK_REAL test_max one = 1.0d0 octant = (CCTK_Equals(domain,"octant")) bitant = (CCTK_Equals(domain,"bitant")) quadrant = (CCTK_Equals(domain,"quadrant")) petsc = (CCTK_Equals(ivp_linear_solver,"petsc")) sor = (CCTK_Equals(ivp_linear_solver,"sor")) jacobi = (CCTK_Equals(ivp_linear_solver,"jacobi")) bamlin = (CCTK_Equals(ivp_linear_solver,"bam")) nonlin = (CCTK_Equals(ivp_linear_solver,"none")) call CCTK_VarIndex (Mlin_index, "ivp::Mlinear") call CCTK_VarIndex (Nsrc_index, "ivp::Nsource") call CCTK_VarIndex (afield_index,"ivp::ivp_psi") if(ivp_linear_bam_precond.eq.1) then write(*,*) 'IVP: calling bam_nonlin as preconditioner to' write(*,*) ' linear solver...' c call ivp_bam_nonlin(GH,my_elliptic_tol) write(*,*) 'BAM Preconditioning NWY' write(*,*) 'Want to put it in? It should be easy ;-)' STOP else ivp_psi = 1.0d0 endif write(*,*)'ivp_linear_ham_tol = ',ivp_linear_ham_tol max_resid = 10.0d0 * ivp_linear_ham_tol icount=0 do while( (max_resid .gt. ivp_linear_ham_tol) .and. . (icount .lt. 10) ) icount=icount+1 ivp_psi_linear_old = ivp_psi Mlinear = ivp_scalarc/8.0d0 - 5.0d0 * ivp_bam_p5 * ivp_psi**4 + . 7.0d0 * ivp_bam_m7 * ivp_psi**(-8) Nsource = 4.0d0 * ivp_bam_p5 * ivp_psi**5 - . 8.0d0 * ivp_bam_m7 * ivp_psi**(-7) #ifdef MAHC write(*,*)'better fix IVP to support MAHC' c$$$ if ((some_hydro .eq. MAHCHYDRO).or. c$$$ . (some_hydro .eq. MAHCHYDRODUST)) then c$$$ if(ivp_matterdens_pow .ne. -8.0d0) then c$$$ write(*,*) 'error:IVP linear_ham_solver: ' c$$$ write(*,*) ' if you want to use ivp_matterdens_pow .ne. -8' c$$$ write(*,*) ' with linear hamiltonian, then you had better' c$$$ write(*,*) ' fix it!' c$$$ STOP c$$$ endif c$$$ Mlin_gf = Mlin_gf + 3.0d0 * ivp_bam_m3 * ivp_psi**(-4) c$$$ Nsrc_gf = Nsrc_gf - 4.0d0 * ivp_bam_m3 * ivp_psi**(-3) c$$$ endif #endif Mlinear = Mlinear + 3.0d0 * ivp_bam_m3 * ivp_psi**(-4) Nsource = Nsource - 4.0d0 * ivp_bam_m3 * ivp_psi**(-3) 1000 continue RelTol(1)=-1 RelTol(2)=-1 RelTol(3)=-1 AbsTol(1)=my_elliptic_tol AbsTol(2)=my_elliptic_tol AbsTol(3)=my_elliptic_tol call CCTK_VarIndex (metpsi_index(1), "einstein::gxx") call CCTK_VarIndex (metpsi_index(2), "einstein::gxy") call CCTK_VarIndex (metpsi_index(3), "einstein::gxz") call CCTK_VarIndex (metpsi_index(4), "einstein::gyy") call CCTK_VarIndex (metpsi_index(5), "einstein::gyz") call CCTK_VarIndex (metpsi_index(6), "einstein::gzz") call CCTK_VarIndex (metpsi_index(7), "einstein::psi") c Linear equations now use a different sign convention. c Malcolm 5/5/2000 Mlinear = -Mlinear Nsource = -Nsource write(*,*)'IVP: calling LinConfMetricSolver' write(*,*)' with solver = ...' ierr = 0.0d0 if(sor.eq.1) then write(*,*)' ...sor' call Ell_LinConfMetricSolver(ierr,cctkGH, . metpsi_index, afield_index, Mlin_index, Nsrc_index, AbsTol, RelTol, . "sor") elseif(petsc.eq.1) then write(*,*)' ...petsc' call Ell_LinConfMetricSolver(ierr,cctkGH, . metpsi_index, afield_index, Mlin_index, Nsrc_index, AbsTol, RelTol, . "petsc") elseif(jacobi.eq.1) then write(*,*)' ...jacobi' call Ell_LinConfMetricSolver(ierr,cctkGH, . metpsi_index, afield_index, Mlin_index, Nsrc_index, AbsTol, RelTol, . "jacobi") elseif(bamlin.eq.1) then write(*,*)' ...BAM linear' write(*,*)'NOTE: Calling BAM through ivp_bam_linear interface' call ivp_bam_linear(cctkGH, my_elliptic_tol) elseif(nonlin.eq.1) then write(*,*)' ...none!' call CCTK_Warn(0,120,"linear_ham_solver","IVP", . "you set ivp_linear_solver = none; choose sor,jacobi,petsc") endif c need to apply inner boundary conditions here bbox_save = cctk_bbox sym(1) = 1 sym(2) = 1 sym(3) = 1 if(octant.eq.1) then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_psi') endif if(quadrant.eq.1)then cctk_bbox(2) = 0 cctk_bbox(4) = 0 cctk_bbox(5) = 0 cctk_bbox(6) = 0 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_psi') endif if(bitant.eq.1)then cctk_bbox(1) = 0 cctk_bbox(2) = 0 cctk_bbox(3) = 0 cctk_bbox(4) = 0 cctk_bbox(6) = 0 call SetCartSymVN(ierr,cctkGH,sym,'ivp::ivp_psi') endif cctk_bbox = bbox_save c calculate residual ivp_psi_linear_old = abs(ivp_psi_linear_old - ivp_psi) c reduce handle call CCTK_ReductionHandle(reduce_handle,"maximum") call CCTK_VarIndex (ivp_psi_linear_old_index, & "ivp::ivp_psi_linear_old") call CCTK_Reduce(ierr, cctkGH, -1, reduce_handle, $ 1,CCTK_VARIABLE_REAL,max_resid, $ 1,ivp_psi_linear_old_index) if (ierr.ne.0) then call CCTK_WARN(1,"Reduction of ivp_psi_linear_old failed"); write(*,*)'ierr = ',ierr c STOP endif write(*,*)'IVP linear_ham_solver:residual ',max_resid enddo c STOP if(icount.ge.100) then write(*,*) 'IVP: Warning: linearization of the ham constraint' write(*,*) ' failed to converge. Must decrease ' write(*,*) ' ivp_linear_ham_tol.' write(*,*) ' Returning with residual: ',max_resid endif return end