Tidying so that the conformal solver isn't used to solve for the conformal

factor, in fact the conformal factor is now only set at the end of
the thorn.

Now if the conformal metric is used, psi will be set to the conformal factor
which was actually solved for. This means that right now that you can't use IDBrillData with conformal derivatives, although it would presumably be
easy enough to numerically calculate them, or to return 1 and 0's instead
just for this case.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDBrillData/trunk@76 a678b1cf-93e1-4b43-a69d-d43939e66649

1 files changed, 45 insertions, 10 deletions

diff --git a/src/finishbrilldata.F b/src/finishbrilldata.F
index 1505d90..9af9fbf 100644
--- a/src/finishbrilldata.F
+++ b/src/finishbrilldata.F
@@ -12,12 +12,13 @@
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
-      subroutine finishbrilldata(CCTK_ARGUMENTS)
+      subroutine IDBrillData_Finish(CCTK_ARGUMENTS)
 
       implicit none
 
       DECLARE_CCTK_ARGUMENTS
       DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
 
       integer i,j,k
       integer nx,ny,nz
@@ -59,7 +60,6 @@ c     Brill metric calculated from q and Psi.
 
                phi = phif(x1,y1)
 
-               psi4 = brillpsi(i,j,k)**4
                e2q  = dexp(2.d0*brillq(rho1,z1,phi))
 
 c              Fudge division by rho^2 on axis. (Physically, y^/rho^2,
@@ -74,20 +74,20 @@ c              The individual coefficients can be read off as
 
                if (rho1.gt.rhofudge) then
 
-                  gxx(i,j,k) = psi4*(e2q + (one - e2q)*y1**2/rho2)
-                  gyy(i,j,k) = psi4*(e2q + (one - e2q)*x1**2/rho2)
-                  gzz(i,j,k) = psi4*e2q
-                  gxy(i,j,k) = - psi4*(one - e2q)*x1*y1/rho2
+                  gxx(i,j,k) = (e2q + (one - e2q)*y1**2/rho2)
+                  gyy(i,j,k) = (e2q + (one - e2q)*x1**2/rho2)
+                  gzz(i,j,k) = e2q
+                  gxy(i,j,k) = - (one - e2q)*x1*y1/rho2
 
                else
 
 c                 This fudge assumes that q = O(rho^2) near the axis. Which
 c                 it should be, or the data will be singular.
 
-                  gxx(i,j,k) = psi4
-                  gyy(i,j,k) = psi4
-                  gzz(i,j,k) = psi4
-                  gxy(i,j,k) = zero
+                  gxx(i,j,k) = 0.0d0
+                  gyy(i,j,k) = 0.0d0
+                  gzz(i,j,k) = 0.0d0
+                  gxy(i,j,k) = 0.0d0
 
                end if
 
@@ -95,6 +95,41 @@ c                 it should be, or the data will be singular.
          end do
       end do
 
+      if (CCTK_EQUALS(metric_type,"static conformal")) then
+
+         conformal_state = 1
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  psi(i,j,k) = brillpsi(i,j,k)
+                  
+               end do
+            end do
+         end do
+
+      else
+
+         conformal_state = 0
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  psi4  = brillpsi(i,j,k)**4
+
+                  gxx(i,j,k) = psi4*gxx(i,j,k)
+                  gyy(i,j,k) = psi4*gyy(i,j,k)
+                  gzz(i,j,k) = psi4*gzz(i,j,k)
+                  gxy(i,j,k) = psi4*gxy(i,j,k)
+                  
+               end do
+            end do
+         end do
+
+      end if
+
 c     In any case,
 
       gxz = zero