Add problem setup for cloud collapse at different metallicities

src/chemistry/Chemistry.hpp

@@ -7,6 +7,7 @@
 //==============================================================================
 /// \file Chemistry.hpp
 /// \brief Defines methods for integrating primordial chemical network using Microphysics
+/// \Authors: Piyush Sharda and Benjamin Wibking
 ///
 
 #include <array>
@@ -19,6 +20,7 @@
 
 #ifdef PRIMORDIAL_CHEM
 #include "actual_eos_data.H"
+#include "actual_network.H"
 #include "burn_type.H"
 #include "eos.H"
 #include "extern_parameters.H"
@@ -39,6 +41,8 @@ template <typename problem_t> auto computeChemistry(amrex::MultiFab &mf, const R
 
 	int num_failed = 0;
 
+	amrex::Real TCMB = 2.73 * (1.0 + network_rp::redshift);
+
 	const BL_PROFILE("Chemistry::computeChemistry()");
 	for (amrex::MFIter iter(mf); iter.isValid(); ++iter) {
 		const amrex::Box &indexRange = iter.validbox();
@@ -90,6 +94,9 @@ template <typename problem_t> auto computeChemistry(amrex::MultiFab &mf, const R
 			// call the EOS to set initial internal energy e
 			eos(eos_input_re, chemstate);
 
+			// set initial Tdust to CMB
+			chemstate.aux[0] = TCMB;
+
 			// do the actual integration
 			// do it in .cpp so that it is not built at compile time for all tests
 			// which would otherwise slow down compilation due to the large RHS file
@@ -114,16 +121,18 @@ template <typename problem_t> auto computeChemistry(amrex::MultiFab &mf, const R
 				insum += inmfracs[nn];
 			}
 
+			// do not normalize ices
 			for (int nn = 0; nn < NumSpec; ++nn) {
-				inmfracs[nn] /= insum;
+
+				if (nn >= network_rp::idx_gas_species) {
+					inmfracs[nn] /= insum;
+				}
 				// update the number densities with conserved mass fractions
 				chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn];
 			}
 
 			// update the number density of electrons due to charge conservation
-			// TODO(psharda): generalize this to other chem networks
-			chemstate.xn[0] = -chemstate.xn[3] - chemstate.xn[7] + chemstate.xn[1] + chemstate.xn[12] + chemstate.xn[6] + chemstate.xn[4] +
-					  chemstate.xn[9] + 2.0 * chemstate.xn[11];
+			balance_charge(chemstate);
 
 			// reconserve mass fractions post charge conservation
 			insum = 0;
@@ -134,7 +143,10 @@ template <typename problem_t> auto computeChemistry(amrex::MultiFab &mf, const R
 			}
 
 			for (int nn = 0; nn < NumSpec; ++nn) {
-				inmfracs[nn] /= insum;
+
+				if (nn >= network_rp::idx_gas_species) {
+					inmfracs[nn] /= insum;
+				}
 				// update the number densities with conserved mass fractions
 				chemstate.xn[nn] = inmfracs[nn] * chemstate.rho / spmasses[nn];
 			}

src/problems/MetalChem/CMakeLists.txt

@@ -0,0 +1,77 @@
+if (AMReX_SPACEDIM EQUAL 3)
+    # Define a custom target that runs the Python script to produce the input perturbations file
+
+    set(microphysics_network_name "metal_chem") #this will override network_name to primordial_chem for this directory only
+    setup_target_for_microphysics_compilation(${microphysics_network_name} "${CMAKE_CURRENT_BINARY_DIR}/")
+
+    #use the BEFORE keyword so that these files get priority in compilation for targets in this directory
+    #this is critical to ensure the correct Microphysics files are linked to metal chem target
+    include_directories(BEFORE ${metal_chem_dirs} "${CMAKE_CURRENT_BINARY_DIR}/" "includes/extern_parameters.H" "includes/network_properties.H")
+
+    add_executable(metalcloud metalcloud.cpp "${QuokkaSourcesNoEOS}" ../../turbulence/TurbDataReader.cpp ../../chemistry/Chemistry.cpp ${metal_chem_sources})
+    target_compile_definitions(metalcloud PUBLIC PRIMORDIAL_CHEM) #this will add #define PRIMORDIAL_CHEM
+
+    if(AMReX_GPU_BACKEND MATCHES "CUDA")
+        setup_target_for_cuda_compilation(metalcloud)
+    endif()
+
+    #need h5py to run perturbations.py file below
+    #purpose of this code is to check whether the h5py python package is installed
+    execute_process(
+        COMMAND Python3::Interpreter -c "h5py"
+        RESULT_VARIABLE EXIT_CODE
+        OUTPUT_QUIET
+    )
+
+    #re-run cmake if this file is changed
+    set_property(DIRECTORY APPEND PROPERTY CMAKE_CONFIGURE_DEPENDS ${CMAKE_SOURCE_DIR}/tests/MetalCloud.in)
+
+    # Read the 'MetalCloud.in' file line by line
+    file(READ ${CMAKE_SOURCE_DIR}/tests/MetalCloud.in pop_in_contents)
+
+    # Split the file contents into lines
+    string(REPLACE "\n" ";" pop_in_lines "${pop_in_contents}")
+
+    # Iterate through the lines and look for lines starting with 'amr.n_cell'
+    foreach(line IN LISTS pop_in_lines)
+        if (line MATCHES "^amr\\.n_cell.*")
+            set(match_line "${line}")
+            break()
+        endif()
+    endforeach()
+
+    if (match_line)
+        message(STATUS "picked line is ${match_line}")
+        # Use a regular expression to find consecutive digit
+        string(REGEX MATCHALL "[0-9]+" digits_list "${match_line}")
+        # Iterate through the matched digits and extract the first one
+        list(GET digits_list 0 first_digits)
+        # Check if matched digits were found
+        if (first_digits)
+           #first_digits give the box size, but you also want box size / 2 as kmax in the script below
+           math(EXPR int_first_digits "${first_digits}")
+           # Check if the conversion was successful
+           if (int_first_digits)
+               # Divide the integer by two
+               math(EXPR result "${int_first_digits} / 2")
+           else()
+               message(FATAL_ERROR "Conversion to integer failed, so cannot find kmax!")
+           endif()
+           message(STATUS "Will create initial perturbations of box size ${first_digits} and kmax ${result}")
+
+        else()
+           message(FATAL_ERROR "No box size found to create initial perturbations for!")
+        endif()
+    else()
+        message(FATAL_ERROR "No matching line found in ${CMAKE_SOURCE_DIR}/tests/MetalCloud.in!")
+    endif()
+
+
+    add_test(NAME ComputePerturbations COMMAND python3 ${CMAKE_SOURCE_DIR}/src/turbulence/perturbation.py --kmin=2 --kmax=${result} --size=${first_digits} --alpha=1.8 --f_solenoidal=0.66667 ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+    add_test(NAME MetalCloud COMMAND metalcloud MetalCloud.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+    set_tests_properties(ComputePerturbations PROPERTIES FIXTURES_SETUP MetalCloud_fixture)
+    set_tests_properties(MetalCloud PROPERTIES FIXTURES_REQUIRED MetalCloud_fixture)
+
+    # AMR test only works on Setonix because Gadi and avatar do not have enough memory per GPU
+    # add_test(NAME MetalCloudAMR COMMAND metalcloud MetalCloud_AMR.in ${QuokkaTestParams} WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}/tests)
+endif()

src/problems/MetalChem/ascent_actions.yaml

@@ -0,0 +1,30 @@
+- 
+  action: "add_pipelines"
+  pipelines: 
+    pl1:
+      f2:
+        type: "slice"
+        params: 
+          point: 
+            x: 0.
+            y: 0.
+            z: 0.
+          normal: 
+            x: 0.0
+            y: 0.0
+            z: 1.0
+- 
+  action: "add_scenes"
+  scenes: 
+    s1: 
+      plots: 
+        p1: 
+          type: "pseudocolor"
+          field: "log_density"
+          pipeline: "pl1"
+      renders:
+        r1:
+          image_prefix: "dens%05d"
+          annotations: "true"
+          camera:
+            zoom: 1.5