[Commits] [svn:einsteintoolkit] GRHydro/trunk/ (Rev. 346)

[rhaas at tapir.caltech.edu]

User: rhaas
Date: 2012/06/05 03:40 PM

Modified:
 /trunk/
  param.ccl, schedule.ccl
 /trunk/src/
  GRHydro_CalcUpdate.F90, GRHydro_ParamCheck.F90

Log:
 GRHydro: Added option to constrain evolution to 1D by
  projecting out non-radial parts of the conserved
  velocity RHSs.
 
 Patch by Christian Reisswig.

File Changes:

Directory: /trunk/src/
======================

File [modified]: GRHydro_CalcUpdate.F90
Delta lines: +64 -0
===================================================================
--- trunk/src/GRHydro_CalcUpdate.F90	2012-06-05 20:32:09 UTC (rev 345)
+++ trunk/src/GRHydro_CalcUpdate.F90	2012-06-05 20:40:27 UTC (rev 346)
@@ -72,6 +72,7 @@
             taurhs(i,j,k) = taurhs(i,j,k) + &
                  (alp_l * tauflux(i-xoffset,j-yoffset,k-zoffset) - &
                  alp_r * tauflux(i,j,k)) * idx
+            
             if(evolve_mhd.ne.0) then
                if(transport_constraints.ne.0) then
                   ! we have to first compute all components of v\crossB = E and
@@ -260,6 +261,7 @@
           taurhs(i,j,k) = taurhs(i,j,k) + &
                (tauflux(i-xoffset,j-yoffset,k-zoffset) - &
                tauflux(i,j,k)) * idx
+          
           if(evolve_mhd.ne.0) then
              Bconsrhs(i,j,k,1) = Bconsrhs(i,j,k,1) + &
                   (Bconsxflux(i-xoffset,j-yoffset,k-zoffset) - &
@@ -319,3 +321,65 @@
 
 end subroutine UpdateCalculation
   
+
+
+
+
+ /*@@
+   @routine    ConstrainSconTo1D
+   @date       Tue  24 14:12 2012
+   @author     Christian Reisswig
+   @desc 
+   Constrains the conserved fluid velocity to radial direction
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @endhistory 
+
+@@*/
+
+
+subroutine ConstrainSconTo1D(CCTK_ARGUMENTS)
+   implicit none
+
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i,j,k
+  CCTK_REAL :: rnorm, rnormI, scon_tmp1, scon_tmp2, scon_tmp3
+
+      !$OMP PARALLEL DO PRIVATE(i,j,k,rnorm,rnormI,scon_tmp1,scon_tmp2,scon_tmp3)
+      do k = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(3) - GRHydro_stencil ! we need to compute Evec on all faces/edges where the fluxes are defined
+        do j = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(2) - GRHydro_stencil
+          do i = GRHydro_stencil + 1 - transport_constraints, cctk_lsh(1) - GRHydro_stencil
+
+            
+            ! Eliminate non-radial fluid velocities to obtain pseudo 1D scheme
+               rnorm = (x(i,j,k)**2 + y(i,j,k)**2 + z(i,j,k)**2)
+               if (rnorm.eq.0) then
+                  rnormI = 0
+               else
+                  rnormI = 1.0d0/rnorm
+               endif
+               
+               scon_tmp1 = (x(i,j,k)*scon(i,j,k,1) &
+                          + y(i,j,k)*scon(i,j,k,2) &
+                          + z(i,j,k)*scon(i,j,k,3)) * rnormI * x(i,j,k)
+               scon_tmp2 = (x(i,j,k)*scon(i,j,k,1) &
+                          + y(i,j,k)*scon(i,j,k,2) &
+                          + z(i,j,k)*scon(i,j,k,3)) * rnormI * y(i,j,k)
+               scon_tmp3 = (x(i,j,k)*scon(i,j,k,1) &
+                          + y(i,j,k)*scon(i,j,k,2) &
+                          + z(i,j,k)*scon(i,j,k,3)) * rnormI * z(i,j,k)
+               
+               scon(i,j,k,1) = scon_tmp1
+               scon(i,j,k,2) = scon_tmp2
+               scon(i,j,k,3) = scon_tmp3
+               
+        end do
+      end do
+    end do
+end subroutine ConstrainSconTo1D
+

File [modified]: GRHydro_ParamCheck.F90
Delta lines: +4 -0
===================================================================
--- trunk/src/GRHydro_ParamCheck.F90	2012-06-05 20:32:09 UTC (rev 345)
+++ trunk/src/GRHydro_ParamCheck.F90	2012-06-05 20:40:27 UTC (rev 346)
@@ -112,6 +112,10 @@
   if (CCTK_Equals(riemann_solver,"HLLE").eq.0.and.coordinates_is_active.ne.0) then
     call CCTK_WARN(0, "There is currently no Riemann solver other than HLLE compatible with multipatch!")
   end if
+  
+  if (constrain_to_1D .eq. 1 .and. coordinates_is_active.eq.1) then
+    call CCTK_WARN(0, "Constrain to 1D option does not work with other than Cartesian coordinates. Ask Christian R. to implement it if required!")
+  endif
 
   if (CCTK_EQUALS(riemann_solver,"Roe").and.CCTK_EQUALS(Bvec_evolution_method,"GRHydro")) then
     call CCTK_PARAMWARN("Roe solver is not implemented yet for MHD")

Directory: /trunk/
==================

File [modified]: param.ccl
Delta lines: +17 -0
===================================================================
--- trunk/param.ccl	2012-06-05 20:32:09 UTC (rev 345)
+++ trunk/param.ccl	2012-06-05 20:40:27 UTC (rev 346)
@@ -612,3 +612,20 @@
 {
 } no
 
+
+
+
+
+############### Testing options ###################################################
+
+
+# This constrains the flux velicity vectors to spherical symmetry.
+# Only the radial contribution to the velocity vectors is considered.
+# This is due to Kotake et al 2012.
+BOOLEAN constrain_to_1D "Set fluid velocities to zero for non-radial motion"
+{
+} no
+
+
+
+

File [modified]: schedule.ccl
Delta lines: +13 -0
===================================================================
--- trunk/schedule.ccl	2012-06-05 20:32:09 UTC (rev 345)
+++ trunk/schedule.ccl	2012-06-05 20:40:27 UTC (rev 346)
@@ -1280,3 +1280,16 @@
     LANG: Fortran
   } "Calculate divB"
 }                           
+
+
+if (constrain_to_1D) {
+  SCHEDULE ConstrainSconTo1D in MoL_PostStepModify
+  {
+    LANG: FORTRAN
+    OPTION: LOCAL
+  } "Constrain conserved fluid velocity to radial direction"
+}
+
+
+
+