From 6150675405532fb739902b920bb2c43431251298 Mon Sep 17 00:00:00 2001 From: miguel Date: Thu, 12 Apr 2001 13:11:38 +0000 Subject: Cleaning up. It is not working correctly yet. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@38 a678b1cf-93e1-4b43-a69d-d43939e66649 --- src/brilldata.F | 95 ++++++++++++++++++++++++++++---------------------- src/setupbrilldata2D.F | 15 ++++---- src/setupbrilldata3D.F | 12 ++++--- 3 files changed, 66 insertions(+), 56 deletions(-) diff --git a/src/brilldata.F b/src/brilldata.F index 4473166..b811ee4 100644 --- a/src/brilldata.F +++ b/src/brilldata.F @@ -4,6 +4,7 @@ #include "cctk_Arguments.h" #include "CactusEinstein/Einstein/src/Einstein.h" +#include "CactusElliptic/EllBase/src/EllBase.h" subroutine brilldata(CCTK_FARGUMENTS) @@ -13,12 +14,28 @@ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer Mlin_index,Nsrc_index,field_index - integer metpsi_index(7) + integer ipsi,iMcoeff,iNcoeff integer ierr CCTK_REAL AbsTol(3),RelTol(3) +c Get indices for grid functions. + + call CCTK_VarIndex(ipsi , "idbrilldata::brillpsi") + if (ipsi .lt. 0) then + call CCTK_WARN(0,"Grid variable index for iphi not found") + end if + + call CCTK_VarIndex(iMcoeff, "idbrilldata::brillMlinear") + if (iMcoeff .lt. 0) then + call CCTK_WARN(0,"Grid variable index for Mcoeff not found") + end if + + call CCTK_VarIndex(iNcoeff, "idbrilldata::brillNsource") + if (iNcoeff .lt. 0) then + call CCTK_WARN(0,"Grid variable index for Ncoeff not found") + end if + c Set up background metric and coefficients for linear solve. if (axisym.eq.1) then @@ -27,64 +44,58 @@ c Set up background metric and coefficients for linear solve. call setupbrilldata3D(CCTK_FARGUMENTS) end if -c Call the Linear Elliptic solver interface to find -c conformal factor. First some preliminaries. - - 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") - - call CCTK_VarIndex(field_index,"IDBrillData::brillpsi") - call CCTK_VarIndex(Mlin_index, "IDBrillData::Mlinear") - call CCTK_VarIndex(Nsrc_index, "IDBrillData::Nsource") +c Tolerances for elliptic solve. AbsTol(1)= brill_thresh - AbsTol(2)= -1 - AbsTol(3)= -1 + AbsTol(2)= brill_thresh + AbsTol(3)= brill_thresh RelTol(1)= -1 RelTol(2)= -1 RelTol(3)= -1 -c Boundaries. +c Elliptic solver. - if (CCTK_EQUALS(brill_bound,"const")) then - call Ell_SetRealKey(ierr,brill_const_v0, - . "EllLinConfMetric::Bnd::Const::V0") - end if + if (axisym.eq.1) then - if (CCTK_EQUALS(brill_bound,"robin")) then - call Ell_SetIntKey(ierr,brill_robin_falloff, - . "EllLinConfMetric::Bnd::Robin::falloff") - call Ell_SetRealKey(ierr,brill_robin_inf, - . "EllLinConfMetric::Bnd::Robin::inf") - end if + call Ell_SetRealKey(ierr, 1.0d0, "EllLinFlat::Bnd::Robin::inf") + call Ell_SetIntKey (ierr, 1, "EllLinFlat::Bnd::Robin::falloff") -c Elliptic solver. + if (CCTK_EQUALS(brill_solver,"sor")) then + call Ell_LinFlatSolver(ierr,cctkGH,ipsi, + . iMcoeff,iNcoeff,AbsTol,RelTol,"sor") + end if - if (CCTK_EQUALS(brill_solver,"sor")) then - call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index, - . field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"sor") - end if + if (CCTK_EQUALS(brill_solver,"petsc")) then + call Ell_LinFlatSolver(ierr,cctkGH,ipsi, + . iMcoeff,iNcoeff,AbsTol,RelTol,"petsc") + end if + + if (CCTK_EQUALS(brill_solver,"bam")) then + call Ell_LinFlatSolver(ierr,cctkGH,ipsi, + . iMcoeff,iNcoeff,AbsTol,RelTol,"bam") + end if + + else - if (CCTK_EQUALS(brill_solver,"petsc")) then - call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index, - . field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"petsc") end if - if (CCTK_EQUALS(brill_solver,"bam")) then - call Ell_LinConfMetricSolver(ierr,cctkGH,metpsi_index, - . field_index,Mlin_index,Nsrc_index,AbsTol,RelTol,"bam") +c Check for errors. + + if (ierr .eq. ELL_SUCCESS) then + call CCTK_INFO("Leaving elliptic solver: solve successful") + else if (ierr .eq. ELL_NOCONVERGENCE) then + call CCTK_INFO("Leaving elliptic solver: solver failed to converge") + else if (ierr .eq. ELL_NOSOLVER) then + call CCTK_INFO("Elliptic solver not found") + else + call CCTK_INFO("Leaving elliptic solver: solve failed") end if c Synchronization and symmetry boundaries. - call CCTK_SyncGroup(cctkGH,"IDBrillData::brillconf") - call CartSymGN(ierr,cctkGH,"IDBrillData::brillconf") + call CCTK_SyncGroup(cctkGH,"idbrilldata::brillconf") + call CartSymGN(ierr,cctkGH,"idbrilldata::brillconf") c Reconstruct physical metric. diff --git a/src/setupbrilldata2D.F b/src/setupbrilldata2D.F index 5643cf0..709370d 100644 --- a/src/setupbrilldata2D.F +++ b/src/setupbrilldata2D.F @@ -1,7 +1,9 @@ + #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" -#include "../../../CactusEinstein/Einstein/src/Einstein.h" + +#include "CactusEinstein/Einstein/src/Einstein.h" subroutine setupbrilldata2D(CCTK_FARGUMENTS) @@ -53,15 +55,10 @@ C Initialize psi. brillpsi = one -C Set up flat metric for the elliptic solve. -C [Delta(g) + Mlinear] psi = Nsource. +C Initialize metric. psi = one - psix = zero - psiy = zero - psiz = zero - gxx = one gyy = one gzz = one @@ -87,7 +84,7 @@ C with a small negative rho, but that should be ok as long as C brillq is even in rho - physically it must be, or the data C will not be regular on the axis. - Mlinear(i,j,k) = - 0.25d0 + brillMlinear(i,j,k) = - 0.25d0 $ *(brillq(rho1,z1+eps,zero) $ + brillq(rho1,z1-eps,zero) $ + brillq(rho1+eps,z1,zero) @@ -100,7 +97,7 @@ C will not be regular on the axis. C Set coefficient Nsource = 0. - Nsource = zero + brillNsource = zero return end diff --git a/src/setupbrilldata3D.F b/src/setupbrilldata3D.F index 22ac6d1..0f6abfa 100644 --- a/src/setupbrilldata3D.F +++ b/src/setupbrilldata3D.F @@ -1,7 +1,9 @@ + #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" -#include "../../../CactusEinstein/Einstein/src/Einstein.h" + +#include "CactusEinstein/Einstein/src/Einstein.h" subroutine setupbrilldata3D(CCTK_FARGUMENTS) @@ -124,7 +126,7 @@ c Find M using centered differences over a small c interval. if (rho1.gt.rhofudge) then - Mlinear(i,j,k) = - 0.25D0/e2q + brillMlinear(i,j,k) = - 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) . + brillq(rho1+eps,z1,phi) @@ -132,7 +134,7 @@ c interval. . - 4.0D0*brillq(rho1,z1,phi)) . / eps**2 else - Mlinear(i,j,k) = - 0.25D0/e2q + brillMlinear(i,j,k) = - 0.25D0/e2q . *(brillq(rho1,z1+eps,phi) . + brillq(rho1,z1-eps,phi) . + two*brillq(eps,z1,phi) @@ -146,7 +148,7 @@ c must always be true otherwise the function c is not regular on the axis. if (rho1.gt.rhofudge) then - Mlinear(i,j,k) = Mlinear(i,j,k) - 0.25D0/rho2 + brillMlinear(i,j,k) = brillMlinear(i,j,k) - 0.25D0/rho2 . *(three*0.25D0*(brillq(rho1,z1,phi+eps) . - brillq(rho1,z1,phi-eps))**2 . + two*(brillq(rho1,z1,phi+eps) @@ -160,7 +162,7 @@ c is not regular on the axis. c Set coefficient Nsource = 0. - Nsource = zero + brillNsource = zero return end -- cgit v1.2.3