FDM Zoom-in simulation

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/Input__DumpTable

@@ -0,0 +1,42 @@
+#Dump ID             Dump Time
+      -6    0.0003124024000000
+      -5    0.0005000000000000
+      -4    0.0010000000000000
+      -3    0.0050000000000000
+      -2    0.0100000000000000
+      -1    0.0500000000000000
+       0    0.0909090920000000
+       1    0.1022642090403268
+       2    0.1142926268126806
+       3    0.1269969647105442
+       4    0.1403806580959482
+       5    0.1544481509848806
+       6    0.1692051241492372
+       7    0.1846587623059424
+       8    0.2008180659019670
+       9    0.2176942139307767
+      10    0.2353009853594673
+      11    0.2536552481440784
+      12    0.2727775264961512
+      13    0.2926926590557780
+      14    0.3134305629162865
+      15    0.3350271209736094
+      16    0.3575252126958853
+      17    0.3809759108571344
+      18    0.4054398686204628
+      19    0.4309889220031104
+      20    0.4577079316383287
+      21    0.4856968849564935
+      22    0.5150732779850334
+      23    0.5459748043551692
+      24    0.5785624241065474
+      25    0.6130240298474579
+      26    0.6495793081783160
+      27    0.6884872598419996
+      28    0.7300595373040358
+      29    0.7746853173110612
+      30    0.8228743661061253
+      31    0.8753084975593694
+      32    0.9327759052183688
+      33    0.9953073615712860
+***************END LINE***************

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/Input__Flag_Interference

@@ -0,0 +1,12 @@
+# Level          QP   Density   PhaseLap   OnlyAtExtrema
+      0        0.03         0        1.0               0
+      1        0.03         0        1.0               0
+      2        0.03         0        1.0               0
+      3        0.03         0        1.0               0
+      4        0.03         0        1.0               0
+      5        0.03         0        1.0               0
+      6        0.03         0        1.0               0
+      7        0.03         0        1.0               0
+      8        0.03         0        1.0               0
+      9        0.03         0        1.0               0
+     10        0.03         0        1.0               0

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/Input__Flag_Rho

@@ -0,0 +1,15 @@
+# Level                         Density
+      0                             8.0
+      1                            32.0
+      2                            64.0
+      3                           256.0
+      4                           330.0
+      5                           430.0
+      6                           512.0
+      7                        262144.0
+      8                       2097152.0
+      9                      16777216.0
+     10                     134217728.0
+     11                    1073741824.0
+     12                    8589934592.0
+     13                   68719476736.0

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/Input__Flag_Spectral

@@ -0,0 +1,13 @@
+# Level   Refinement   Derefinement (derefinement currently not functional)
+      0          1e0             -1
+      1          1e0             -1
+      2          1e0             -1
+      3          1e0             -1
+      4          1e0             -1
+      5          1e0             -1
+      6          1e0             -1
+      7          1e0             -1
+      8          1e0             -1
+      9          1e0             -1
+     10          1e0             -1
+     11          1e0             -1

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/Input__TestProb

@@ -0,0 +1,3 @@
+# problem-specific runtime parameters
+LSS_InitMode         2       # initialization mode: (1=density-only with constant phase, 2=real and imaginary parts or density and phase) [1]
+ZoomIn_MaxLvOutside  3       # if (OPT__REIGON_FLAG), maximum refinement level outside of the zoom-in box. Required: ELBDM_FIRST_WAVE_LEVEL <= ZoomIn_MaxLvOutside

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/Input__ZoominBox

@@ -0,0 +1,8 @@
+#ScaleFactor    LX(UNIT_L)    LY(UNIT_L)    LZ(UNIT_L)   CenterX(UNIT_L)   CenterY(UNIT_L)   CenterZ(UNIT_L)
+        0.77          1.18          1.13          1.26           0.29491           9.52254           8.27318
+	0.58          1.80          1.66          1.94           0.29491           9.52254           8.27318
+	0.43          2.30          2.00          2.46           0.29491           9.52254           8.27318
+	0.31          2.70          2.31          2.87           0.29491           9.52254           8.27318
+	0.22          3.20          2.53          3.22           0.29491           9.52254           8.27318
+	0.14          3.60          2.70          3.50           0.29491           9.52254           8.27318
+	0.05          4.00          2.90          3.89           0.29491           9.52254           8.27318

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/README

@@ -0,0 +1,50 @@
+Compilation flags:
+========================================
+Enable : MODEL=ELBDM, ELBDM_SCHEME=ELBDM_HYBRID, GRAVITY, COMOVING
+Disable: PARTICLE
+
+
+Default setup:
+========================================
+1. ELBDM_MASS              2.0e-24 (eV/c^2)
+   A_INIT                  9.900990099e-3          # initial scale factor
+   OMEGA_M0                0.3158230904284232      # omega matter at the present time
+   HUBBLE0                 0.6732117               # dimensionless Hubble parameter (currently only for converting ELBDM_MASS to code units)
+   BOX_SIZE                10. (Mpc/h)
+
+2. MAX_LEVEL               7
+   OPT__FLAG_SPECTRAL      1
+   OPT__FLAG_INTERFERENCE  1
+   --> Input__Flag_Rho          = example/input/Input__Flag_Rho
+       Input__Flag_Spectral     = example/input/Input__Flag_Spectral
+       Input__Flag_Interference = example/input/Input__Flag_Interference
+
+3. OPT__FLAG_REGION        1			   # Enable zoom-in simulation
+
+
+Note:
+========================================
+1. Cosmological ZOOM-IN large-scale structure simulations using hybrid scheme
+
+2. Quickstart:
+   1. Download the IC file "sh download_IC_lighthalo.sh" for testing.
+   2. Compile GAMER with the "Makefile" generated by "generate_make.sh" and run simulation to redshift 0 with default settings.
+   3. Total data size is 74G. If less snapshots are preferred for small storage space, please remove some rows in Input__DumpTable
+
+3. details on Zoom-in Simulation
+   1. For details on running a hybrid simulation, please see "example/test_problem/ELBDM/LSS_Hybrid/README"
+   2. Input__ZoomInBox controls the volume(LX, LY, LZ) and center(CenterX, CenterY, CenterZ) of the time dependent zoom-in box.
+      Note that we read the scale factor in descending order.
+      The UNIT_L in cosmological simulation refers to cMpc/h.
+   3. Input__TestProb controls the maximum level outside of the volume (ZoomIn_MaxLvOutside).
+   4. It is required to set ELBDM_FIRST_WAVE_LEVEL <= ZoomIn_MaxLvOutside, to avoid problems of the fluid-wave interface at the zoom-in boundary (Chan+24, under review)
+   5. OPT__FLAG_REGION must be enabled.
+
+4. Some yt visualization scripts are put in "plot_script"
+
+5. Simulation results in (Chan+24, submitted to MNRAS)
+   1. "download_IC_heavyhalo.sh" downloads the IC file
+   2. NX0_TOT_X   256
+      NX0_TOT_Y   256
+      NX0_TOT_Z   256
+      ELBDM_MASS  2.0e-23

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/clean.sh

@@ -0,0 +1,5 @@
+rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Dump Record__MemInfo Record__L1Err \
+      Record__Conservation Data* stderr stdout log XYslice* YZslice* XZslice* Xline* Yline* Zline* \
+      Diag* BaseXYslice* BaseYZslice* BaseXZslice* BaseXline* BaseYline* BaseZline* BaseDiag* \
+      PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \
+      Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance Record__Hybrid

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/download_IC_heavyhalo.sh

@@ -0,0 +1,13 @@
+#!/bin/bash
+
+API_URL="https://girder.hub.yt/api/v1"
+FILE_ID="66ff5e05999605c485c8d67a"
+LOCAL_FILE="Zoomin_IC"
+
+# download
+girder-cli --api-url ${API_URL} download --parent-type item ${FILE_ID} ${LOCAL_FILE}
+
+# unzip
+tar zxvf ${LOCAL_FILE}/Zoomin_IC_heavyhalo.tar.gz
+rm -r ${LOCAL_FILE}
+ln -fs UM_IC_hybrid_N1024 UM_IC

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/download_IC_lighthalo.sh

@@ -0,0 +1,13 @@
+#!/bin/bash
+
+API_URL="https://girder.hub.yt/api/v1"
+FILE_ID="66ff5e15999605c485c8d67d"
+LOCAL_FILE="Zoomin_IC"
+
+# download
+girder-cli --api-url ${API_URL} download --parent-type item ${FILE_ID} ${LOCAL_FILE}
+
+# unzip
+tar zxvf ${LOCAL_FILE}/Zoomin_IC_lighthalo.tar.gz
+rm -r ${LOCAL_FILE}
+ln -fs UM_IC_hybrid_N256_third UM_IC

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/generate_make.sh

@@ -0,0 +1,7 @@
+# This script should run in the same directory as configure.py
+
+PYTHON=python3
+
+${PYTHON} configure.py --machine=eureka_intel --mpi=true --hdf5=true --fftw=FFTW3 --gpu=true \
+                       --model=ELBDM --elbdm_scheme=ELBDM_HYBRID --wave_scheme=WAVE_GRAMFE --gramfe_scheme=GRAMFE_MATMUL \
+                       --gravity=true --comoving=true --gsl=true --spectral_interpolation=true

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/plot_script/Compute_profiles.py

@@ -0,0 +1,60 @@
+#!/usr/bin/env python3.9
+################################################################################
+
+# This code outputs the density, cumulative mass, circular velocity and velocity dispersion profiles
+# of a target halo. The code uses an initial estimate (center_first_guess) to search the vicinity within
+# a specified radius (vicinity) for the coordinates of the maximum density,
+# corresponding to the center of the target halo.
+
+# carefully adjust vicinity for performance.
+
+# Example usage after running simulation with low resolution IC:
+#    python3 Compute_profiles.py -s 36 -e 36
+
+# Outputs:
+#   Halo_Parameter
+#   prof_dens/Data_xxxxxx_0_profile_data
+#   prof_mass/Data_xxxxxx_0_mass_accumule
+#   prof_dens/Data_xxxxxx_0_profile_data
+#   prof_circular_vel/Data_xxxxxx_0_circular_velocity
+#   prof_veldisp/Data_xxxxxx_0_veldisp_haloRestFrame
+
+################################################################################
+import argparse
+import numpy as np
+import yt
+import sys
+from Profile_Functions import *
+
+
+# load the command-line parameters
+parser = argparse.ArgumentParser( description='Plot profile and output Halo_parameter' )
+
+parser.add_argument( '-s', action='store', required=True,  type=int, dest='idx_start',
+                      help='first data index' )
+parser.add_argument( '-e', action='store', required=True,  type=int, dest='idx_end',
+                      help='last data index' )
+parser.add_argument( '-d', action='store', required=False, type=int, dest='didx',
+                      help='delta data index [%(default)d]', default=1 )
+
+# take note
+print( '\nCommand-line arguments:' )
+print( '-------------------------------------------------------------------' )
+print( ' '.join(map(str, sys.argv)) )
+print( '-------------------------------------------------------------------\n' )
+
+args        = parser.parse_args()
+idx_start   = args.idx_start
+idx_end     = args.idx_end
+didx        = args.didx
+
+halo_id      = 0   # pick an id for your target halo
+
+for file_id in range( idx_start, idx_end+1, didx ):
+
+    ds                 = yt.load( "../Data_%06d"%file_id )
+    center_first_guess = np.array( [0.295, 9.522, 8.27] ) # in cMpc/h. First guess for target halo of low resolution IC at z=0
+    vicinity           = 0.3                              # radius in cMpc/h
+    compute_profile( ds, center_first_guess, vicinity, halo_id, '.' )
+
+print( "Done !" )

example/test_problem/ELBDM/LSS_Hybrid_Zoomin/plot_script/ELBDM_DerivedField.py

@@ -0,0 +1,127 @@
+import yt
+import numpy as np
+
+# define the derived fields
+###########################
+# psi = R + iI = f*e^(iS)
+
+# amplitude
+# = sqrt(rho) = sqrt(R^2 + I^2)
+def _f( field, data ):
+   return data["Dens"]**0.5
+
+# real part
+# = sqrt(rho)*cos(S) = R
+def _Real( field, data ):
+   return data["f"]*np.cos( data["Phase"] )
+
+# imaginary part
+# = sqrt(rho)*sin(S) = I
+def _Imag( field, data ):
+   return data["f"]*np.sin( data["Phase"] )
+
+# phase
+# = arctan(I/R)
+def _S( field, data ):
+   return np.arctan2( data["Imag"], data["Real"] )
+
+def _phasemod2( field, data ):
+   return np.arctan2( np.sin(data["Phase"]), np.cos(data["Phase"]) )
+
+# momentum density x-component
+# = rho*v_x = (R*dI/dx - I*dR/dx)*hbar/m
+def _j_x( field, data ):
+   ELBDM_ETA = data.ds.parameters['ELBDM_Mass']*data.ds.units.code_mass/data.ds.units.reduced_planck_constant
+   return (data["Real"]*data["Imag_gradient_x"] - data["Imag"]*data["Real_gradient_x"])/ELBDM_ETA
+
+# momentum density y-component
+# = rho*v_y = (R*dI/dy - I*dR/dy)*hbar/m
+def _j_y( field, data ):
+   ELBDM_ETA = data.ds.parameters['ELBDM_Mass']*data.ds.units.code_mass/data.ds.units.reduced_planck_constant
+   return (data["Real"]*data["Imag_gradient_y"] - data["Imag"]*data["Real_gradient_y"])/ELBDM_ETA
+
+# momentum density z-component
+# = rho*v_z = (R*dI/dz - I*dR/dz)*hbar/m
+def _j_z( field, data ):
+   ELBDM_ETA = data.ds.parameters['ELBDM_Mass']*data.ds.units.code_mass/data.ds.units.reduced_planck_constant
+   return (data["Real"]*data["Imag_gradient_z"] - data["Imag"]*data["Real_gradient_z"])/ELBDM_ETA
+
+# momentum density magnitude
+# = rho*|v| = |R*grad(I) - I*grad(R)|*hbar/m
+def _j( field, data ):
+   return (data["j_x"]**2 + data["j_y"]**2 + data["j_z"]**2)**0.5
+
+# bulk velocity x-component
+# = (dS/dx)*hbar/m = (R*dI/dx - I*dR/dx)/rho*hbar/m
+def _v_x( field, data ):
+   return data["j_x"]/data["Dens"]
+
+# bulk velocity y-component
+# = (dS/dy)*hbar/m = (R*dI/dy - I*dR/dy)/rho*hbar/m
+def _v_y( field, data ):
+   return data["j_y"]/data["Dens"]
+
+# bulk velocity z-component
+# = (dS/dz)*hbar/m = (R*dI/dz - I*dR/dz)/rho*hbar/m
+def _v_z( field, data ):
+   return data["j_z"]/data["Dens"]
+
+# bulk velocity magnitude
+# = |grad(S)|*hbar/m = |(R*grad(I) - I*grad(R))|/rho*hbar/m
+def _v( field, data ):
+   return (data["v_x"]**2 + data["v_y"]**2 + data["v_z"]**2)**0.5
+
+###########################
+
+
+def Add_derived_fields( ds ):
+
+   ds.add_field( ("gamer","f"),
+                 function=_f, units="code_mass**0.5/code_length**1.5",
+                 display_name=r"$f$", sampling_type="cell" )
+   ds.add_field( ("gamer","Real"),
+                 function=_Real, units="code_mass**0.5/code_length**1.5",
+                 display_name=r"$Real$", sampling_type="cell" )
+   ds.add_field( ("gamer","Imag"),
+                 function=_Imag, units="code_mass**0.5/code_length**1.5",
+                 display_name=r"$Imag$", sampling_type="cell" )
+   ds.add_field( ("gamer","S"),
+                 function=_S, units="dimensionless",
+                 display_name=r"$S$", sampling_type="cell" )
+   ds.add_field( ("gamer", "PhaseMod2"),
+                 function=_phasemod2,
+                 sampling_type="local", units="" )
+
+   # gradient fields
+   Grad_R = ds.add_gradient_fields( ("gamer","Real") )
+   Grad_I = ds.add_gradient_fields( ("gamer","Imag") )
+   Grad_f = ds.add_gradient_fields( ("gamer","f") )
+   Grad_S = ds.add_gradient_fields( ("gamer","S") )
+
+   # momentum density
+   ds.add_field( ("gamer","j_x"),
+                 function=_j_x, units="code_mass/(code_length**2*code_time)",
+                 display_name=r"$j_x$", sampling_type="cell" )
+   ds.add_field( ("gamer","j_y"),
+                 function=_j_y, units="code_mass/(code_length**2*code_time)",
+                 display_name=r"$j_y$", sampling_type="cell" )
+   ds.add_field( ("gamer","j_z"),
+                 function=_j_z, units="code_mass/(code_length**2*code_time)",
+                 display_name=r"$j_z$", sampling_type="cell" )
+   ds.add_field( ("gamer","j"),
+                 function=_j, units="code_mass/(code_length**2*code_time)",
+                 display_name=r"$|j|$", sampling_type="cell" )
+
+   # velocity
+   ds.add_field( ("gamer","v_x"),
+                 function=_v_x, units="code_length/code_time",
+                 display_name=r"$v_x$", sampling_type="cell" )
+   ds.add_field( ("gamer","v_y"),
+                 function=_v_y, units="code_length/code_time",
+                 display_name=r"$v_y$", sampling_type="cell" )
+   ds.add_field( ("gamer","v_z"),
+                 function=_v_z, units="code_length/code_time",
+                 display_name=r"$v_z$", sampling_type="cell" )
+   ds.add_field( ("gamer","v"),
+                 function=_v, units="code_length/code_time",
+                 display_name=r"$|v|$", sampling_type="cell" )