MHD (#307)

This is a first implementation of a contrained transport, corner transport upwind MHD solver based on Miniati & Colella 2011.  At the moment, it is not tested with AMR.

Exec/Make.Castro

@@ -148,7 +148,16 @@ ifeq ($(USE_ALL_CASTRO), TRUE)
            Source/scf Util/model_parser Util/conservative_interpolate
 endif
 
+
 # add / define any special physics we need
+
+ifeq ($(USE_MHD), TRUE)
+  Bdirs += Source/mhd
+  DEFINES += -DMHD
+else
+  Bdirs += Source/hydro
+endif
+
 ifeq ($(USE_GRAV), TRUE)
   Bdirs += Source/gravity
   DEFINES += -DGRAVITY

Exec/gravity_tests/DustCollapse/GNUmakefile

@@ -13,6 +13,8 @@ USE_MPI    = TRUE
 
 USE_GRAV   = TRUE
 
+USE_MHD = FALSE
+
 USE_PROB_PARAMS = TRUE
 
 GPU_COMPATIBLE_PROBLEM = TRUE

Exec/gravity_tests/DustCollapse/Prob_nd.F90

@@ -188,3 +188,57 @@ subroutine ca_initdata(lo, hi, &
   end subroutine ca_initdata
 
 end module initdata_module
+
+#ifdef MHD
+subroutine ca_initmag(level, time, lo, hi, &
+                      nbx, mag_x, bx_lo, bx_hi, &
+                      nby, mag_y, by_lo, by_hi, &
+                      nbz, mag_z, bz_lo, bz_hi, &
+                      delta, xlo, xhi)
+
+  use probdata_module, only : B_x, B_y, B_z
+  use prob_params_module
+  use amrex_fort_module, only : rt => amrex_real
+
+  implicit none
+
+  integer :: level, nbx, nby, nbz
+  integer :: lo(3), hi(3)
+  integer :: bx_lo(3), bx_hi(3)
+  integer :: by_lo(3), by_hi(3)
+  integer :: bz_lo(3), bz_hi(3)
+  real(rt) :: xlo(3), xhi(3), time, delta(3)
+
+  real(rt) :: mag_x(bx_lo(1):bx_hi(1), bx_lo(2):bx_hi(2), bx_lo(3):bx_hi(3), nbx)
+  real(rt) :: mag_y(by_lo(1):by_hi(1), by_lo(2):by_hi(2), by_lo(3):by_hi(3), nby)
+  real(rt) :: mag_z(bz_lo(1):bz_hi(1), bz_lo(2):bz_hi(2), bz_lo(3):bz_hi(3), nbz)
+
+  real(rt) :: xcen, ycen, zcen
+  integer  :: i, j, k
+
+
+  do k = lo(3), hi(3)
+     do j = lo(2), hi(2)
+        do i = lo(1), hi(1)+1
+              mag_x(i,j,k,1) = B_x
+        enddo
+     enddo
+  enddo
+  do k = lo(3), hi(3)
+     do j = lo(2), hi(2)+1
+        do i = lo(1), hi(1)
+           mag_y(i,j,k,1) = B_y
+         enddo
+     enddo
+  enddo
+
+  do k = lo(3), hi(3)+1
+     do j = lo(2), hi(2)
+        do i = lo(1), hi(1)
+           mag_z(i,j,k,1) = B_z
+        enddo
+     enddo
+  enddo
+
+end subroutine ca_initmag
+#endif

Exec/gravity_tests/DustCollapse/_prob_params

@@ -28,4 +28,9 @@ center_y               real            0.0_rt         y
 
 center_z               real            0.0_rt         y
 
-nsub                   integer         5              y
+nsub                   integer         5              y
+
+B_x                    real            0.0e0_rt       y
+B_y                    real            0.0e0_rt       y
+B_z                    real            0.0e0_rt       y
+

Exec/gravity_tests/DustCollapse/bc_fill_nd.F90

@@ -285,4 +285,71 @@ subroutine denfill(lo, hi, adv, adv_lo, adv_hi, domlo, domhi, delta, xlo, time,
 
   end subroutine denfill
 
+#ifdef MHD
+  subroutine face_fillx(lo, hi, var, var_lo, var_hi, domlo, domhi, delta, xlo, time, bc) bind(C, name="face_fillx")
+
+    use amrex_fort_module, only : rt => amrex_real
+    use fc_fill_module
+
+    implicit none
+
+    integer,  intent(in   ) :: lo(3), hi(3)
+    integer,  intent(in   ) :: var_lo(3), var_hi(3)
+    integer,  intent(in   ) :: domlo(3), domhi(3)
+    integer,  intent(in   ) :: bc(AMREX_SPACEDIM,2)
+    real(rt), intent(in   ) :: delta(3), xlo(3)
+    real(rt), intent(inout) :: var(var_lo(1):var_hi(1), var_lo(2):var_hi(2), var_lo(3):var_hi(3))
+    real(rt), intent(in   ), value :: time
+    integer dir
+
+    dir = 1
+
+    call filfc(var,var_lo(1),var_lo(2),var_lo(3),var_hi(1),var_hi(2),var_hi(3),domlo,domhi,delta,xlo,bc,dir)
+
+  end subroutine face_fillx
+
+
+  subroutine face_filly(lo, hi, var, var_lo, var_hi, domlo, domhi, delta, xlo, time, bc) bind(C, name="face_filly")
+
+    use amrex_fort_module, only : rt => amrex_real
+    use fc_fill_module
+
+    implicit none
+    integer,  intent(in   ) :: lo(3), hi(3)
+    integer,  intent(in   ) :: var_lo(3), var_hi(3)
+    integer,  intent(in   ) :: domlo(3), domhi(3)
+    integer,  intent(in   ) :: bc(AMREX_SPACEDIM,2)
+    real(rt), intent(in   ) :: delta(3), xlo(3)
+    real(rt), intent(inout) :: var(var_lo(1):var_hi(1), var_lo(2):var_hi(2), var_lo(3):var_hi(3))
+    real(rt), intent(in   ), value :: time
+    integer dir
+
+    dir = 2
+
+    call filfc(var,var_lo(1),var_lo(2),var_lo(3),var_hi(1),var_hi(2),var_hi(3),domlo,domhi,delta,xlo,bc,dir)
+
+  end subroutine face_filly
+
+  subroutine face_fillz(lo, hi, var, var_lo, var_hi, domlo, domhi, delta, xlo, time, bc) bind(C, name="face_fillz")
+
+    use amrex_fort_module, only : rt => amrex_real
+    use fc_fill_module
+
+    implicit none
+
+    integer,  intent(in   ) :: lo(3), hi(3)
+    integer,  intent(in   ) :: var_lo(3), var_hi(3)
+    integer,  intent(in   ) :: domlo(3), domhi(3)
+    integer,  intent(in   ) :: bc(AMREX_SPACEDIM,2)
+    real(rt), intent(in   ) :: delta(3), xlo(3)
+    real(rt), intent(inout) :: var(var_lo(1):var_hi(1), var_lo(2):var_hi(2), var_lo(3):var_hi(3))
+    real(rt), intent(in   ), value :: time
+    integer dir
+    dir = 3
+
+    call filfc(var,var_lo(1),var_lo(2),var_lo(3),var_hi(1),var_hi(2),var_hi(3),domlo,domhi,delta,xlo,bc,dir)
+
+  end subroutine face_fillz
+#endif
+
 end module bc_fill_module

Exec/gravity_tests/DustCollapse/inputs_3d_monopole_mhd_test

@@ -0,0 +1,68 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 10
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0 0 0
+geometry.coord_sys   = 0
+geometry.prob_lo     = 0.    0.    0.
+geometry.prob_hi     = 7.5e8 7.5e8 7.5e8
+
+amr.n_cell           =  64    64    64
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  3 3 3
+castro.hi_bc       =  2 2 2
+
+# WHICH PHYSICS
+castro.do_hydro      = 1
+castro.do_react      = 0
+castro.add_ext_src   = 0
+castro.do_sponge     = 1
+castro.ppm_type      = 1
+
+# MHD
+castro.mhd_plm_slope = 1
+castro.use_flattening = 0
+#castro.castro.apply_sources_consecutively = 1
+
+castro.do_grav       = 1
+gravity.gravity_type = MonopoleGrav
+gravity.drdxfac      = 4
+
+# RELEVANT FOR EOS
+castro.small_temp    = 1.e-3
+castro.small_dens    = 1.e-6
+
+# TIME STEP CONTROL
+castro.cfl            = 0.5     # cfl number for hyperbolic system
+castro.init_shrink    = 0.1     # scale back initial timestep
+castro.change_max     = 1.05    # scale back initial timestep
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 4       # block factor in grid generation
+amr.max_grid_size   = 32
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = chk_3d_  # root name of checkpoint file
+amr.check_int       = 100000        # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = plt_3d_  # root name of plotfile
+amr.plot_int        = 100000        # number of timesteps between plotfiles
+amr.derive_plot_vars = NONE
+
+#PROBIN FILENAME
+amr.probin_file = probin.octant

Exec/gravity_tests/DustCollapse/inputs_3d_poisson_mhd_test

@@ -0,0 +1,72 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 8
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic = 0 0 0
+geometry.coord_sys   = 0
+geometry.prob_lo     = 0.    0.    0.
+geometry.prob_hi     = 1.5e9 1.5e9 1.5e9
+
+amr.n_cell           =  64  64  64
+
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+# 0 = Interior           3 = Symmetry
+# 1 = Inflow             4 = SlipWall
+# 2 = Outflow            5 = NoSlipWall
+# >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
+castro.lo_bc       =  2 2 2
+castro.hi_bc       =  2 2 2
+
+# WHICH PHYSICS
+castro.do_hydro      = 1
+castro.do_react      = 0
+castro.add_ext_src   = 0
+castro.do_sponge     = 1
+castro.ppm_type      = 1
+castro.do_grav       = 1
+gravity.gravity_type = PoissonGrav
+gravity.drdxfac      = 4
+
+#gravity.v = 1
+#mg.v = 2
+
+# RELEVANT FOR EOS
+castro.small_temp    = 1.e-3
+castro.small_dens    = 1.e-6
+
+# MHD 
+castro.mhd_plm_slope = 1
+castro.use_flattening = 0
+#castro.castro.apply_sources_consecutively = 1
+
+
+# TIME STEP CONTROL
+castro.cfl            = 0.5     # cfl number for hyperbolic system
+castro.init_shrink    = 0.03    # scale back initial timestep
+castro.change_max     = 1.05    # scale back initial timestep
+
+# DIAGNOSTICS & VERBOSITY
+castro.sum_interval   = 1       # timesteps between computing mass
+castro.v              = 1       # verbosity in Castro.cpp
+amr.v                 = 1       # verbosity in Amr.cpp
+#amr.grid_log        = grdlog  # name of grid logging file
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 0       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 4       # block factor in grid generation
+amr.max_grid_size   = 32
+amr.n_error_buf     = 2 2 2 2 # number of buffer cells in error est
+
+# CHECKPOINT FILES
+amr.check_file      = chk_3d_  # root name of checkpoint file
+amr.check_int       = 8        # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_file       = plt_3d_  # root name of plotfile
+amr.plot_int        = 8        # number of timesteps between plotfiles
+amr.derive_plot_vars = NONE
+
+#PROBIN FILENAME
+amr.probin_file = probin_3d_poisson_regtest

Exec/hydro_tests/Sedov/GNUmakefile

@@ -10,6 +10,8 @@ COMP	   = gnu
 USE_MPI    = TRUE
 USE_OMP    = FALSE
 
+USE_MHD    = FALSE
+
 GPU_COMPATIBLE_PROBLEM = TRUE
 USE_PROB_PARAMS = TRUE
 

Exec/hydro_tests/Sedov/Prob_nd.F90

@@ -276,3 +276,60 @@ subroutine ca_initdata(lo, hi, &
   enddo
 
 end subroutine ca_initdata
+
+#ifdef MHD
+subroutine ca_initmag(level, time, lo, hi, &
+                      nbx, mag_x, bx_lo, bx_hi, &
+                      nby, mag_y, by_lo, by_hi, &
+                      nbz, mag_z, bz_lo, bz_hi, &
+                      delta, xlo, xhi)
+
+  use probdata_module
+  use prob_params_module, only : center
+  use amrex_fort_module, only : rt => amrex_real
+
+  implicit none
+
+  integer :: level, nbx, nby, nbz
+  integer :: lo(3), hi(3)
+  integer :: bx_lo(3), bx_hi(3)
+  integer :: by_lo(3), by_hi(3)
+  integer :: bz_lo(3), bz_hi(3)
+  real(rt) :: xlo(3), xhi(3), time, delta(3)
+
+  real(rt) :: mag_x(bx_lo(1):bx_hi(1), bx_lo(2):bx_hi(2), bx_lo(3):bx_hi(3), nbx)
+  real(rt) :: mag_y(by_lo(1):by_hi(1), by_lo(2):by_hi(2), by_lo(3):by_hi(3), nby)
+  real(rt) :: mag_z(bz_lo(1):bz_hi(1), bz_lo(2):bz_hi(2), bz_lo(3):bz_hi(3), nbz)
+
+  real(rt) :: xcen, ycen, zcen
+  integer  :: i, j, k
+
+
+  do k = lo(3), hi(3)
+     do j = lo(2), hi(2)
+        do i = lo(1), hi(1)+1
+           mag_x(i,j,k,1) = B_x
+        enddo
+     enddo
+  enddo
+
+  do k = lo(3), hi(3)
+     do j = lo(2), hi(2)+1
+        do i = lo(1), hi(1)
+           mag_y(i,j,k,1) = B_y
+        enddo
+     enddo
+  enddo
+
+  do k = lo(3), hi(3)+1
+     do j = lo(2), hi(2)
+        do i = lo(1), hi(1)
+           mag_z(i,j,k,1) = B_z
+        enddo
+     enddo
+  enddo
+
+end subroutine ca_initmag
+
+
+#endif

Exec/hydro_tests/Sedov/_prob_params

@@ -19,3 +19,8 @@ nsub                   integer              4                            y
 
 e_exp                  real                 0.0_rt
 
+
+B_x                    real                 0.0_rt                       y
+B_y                    real                 0.0_rt                       y
+B_z                    real                 0.0_rt                       y
+