diff --git a/data/cav5_comb4/3D/scalability/mkmsh b/data/cav5_comb4/3D/scalability/mkmsh index f9e7e0a..68c1258 100755 --- a/data/cav5_comb4/3D/scalability/mkmsh +++ b/data/cav5_comb4/3D/scalability/mkmsh @@ -18,6 +18,7 @@ SAMPFAC="4" OUTPUT="actii" MKLINK="false" OVRWT="false" +NELEM="0" while [[ $# -gt 0 ]]; do arg=$1 @@ -65,6 +66,9 @@ while [[ $# -gt 0 ]]; do --output=*) OUTPUT=${arg#*=} ;; + --nelem=*) + NELEM=${arg#*=} + ;; --link) MKLINK="true" ;; @@ -94,7 +98,9 @@ printf "Cavity factor: ${CAVFAC}\n" printf "Injector factor: ${INJFAC}\n" printf "Shear layer factor: ${SHEARFAC}\n" printf "Sample factor: ${SAMPFAC}\n" + if [[ "${OVRWT}" == "true" ]]; then + rm -f ${OUTPUT}_${NELEM}.msh rm -f ${OUTPUT}.msh MKLINK="true" fi @@ -106,6 +112,13 @@ if [[ "${MKLINK}" == "true" ]]; then fi printf "Make link ${OUTPUT}.msh for finished mesh.\n" fi +if [[ ! -f ${OUTPUT}_${NELEM}.msh ]]; then + date + set -x + rm -f tmp.msh + gmsh -setnumber size ${SIZE} -setnumber blratio ${BLRAT} -setnumber cavityfac ${CAVFAC} -setnumber isofac ${ISOFAC} -setnumber samplefac ${SAMPFAC} -setnumber injectorfac ${INJFAC} -setnumber blratiocavity ${BLCAVRAT} -setnumber blratioinjector ${BLINJRAT} -setnumber blratiosample ${BLSAMPRAT} -setnumber blratiosurround ${BLSURRRAT} -setnumber shearfac ${SHEARFAC} -o tmp.msh -nopopup -format msh2 ./actii_from_brep.geo -3 -nt ${NCPUS} + set +x + date date set -x @@ -120,13 +133,11 @@ if [ -f "tmp.msh" ]; then printf "Mesh creation failed.\n" exit 1 fi - mv tmp.msh ${OUTPUT}_${NCELLS}.msh - printf "Created ${OUTPUT}_${NCELLS}.msh with ${NCELLS} tets.\n" +else + printf "Expected mesh file already exists, skipping mesh creation.\n" if [[ "${MKLINK}" == "true" ]]; then - ln -s ${OUTPUT}_${NCELLS}.msh ${OUTPUT}.msh + ln -s ${OUTPUT}_${NELEM}.msh ${OUTPUT}.msh fi -else - printf "Mesh creation failed." - exit 1 fi + diff --git a/data/cav5_comb4/3D/scalability/sizes.txt b/data/cav5_comb4/3D/scalability/sizes.txt index 2ea556b..411fbcb 100644 --- a/data/cav5_comb4/3D/scalability/sizes.txt +++ b/data/cav5_comb4/3D/scalability/sizes.txt @@ -9,4 +9,4 @@ | 7.45 | 1,536k | 1,536,571 | | 5.8545 | 3,072k | 3,071,835 | | 4.61 | 6,144k | 6,142,150 | -| 3.6389 | 12,288k | 12,283,032 | +| 3.6389 | 12,288k | 12,283,032 | \ No newline at end of file diff --git a/flash_test/driver.py b/flash_test/driver.py new file mode 120000 index 0000000..3f406d1 --- /dev/null +++ b/flash_test/driver.py @@ -0,0 +1 @@ +../driver.py \ No newline at end of file diff --git a/flash_test/flash_test.bsub.sh b/flash_test/flash_test.bsub.sh new file mode 100755 index 0000000..bca2e1a --- /dev/null +++ b/flash_test/flash_test.bsub.sh @@ -0,0 +1,20 @@ +#! /bin/bash +#BSUB -nnodes 1 +#BSUB -G uiuc +#BSUB -W 360 +#BSUB -J flash9 +#BSUB -q pbatch +#BSUB -o flash9.out + +source ../emirge/config/activate_env.sh +source ../emirge/mirgecom/scripts/mirge-testing-env.sh + +set -x +$MIRGE_MPI_EXEC -n 1 $MIRGE_PARALLEL_SPAWNER python -u -O -m mpi4py driver.py -c ${casename} -i run_params.yaml -r restart_data/prediction-000965000 --log --lazy +set +x + +# flash5 000483000 +# flash6 000604000 +# flash7 000726000 +# flash8 000847000 +# flash9 000965000 diff --git a/flash_test/run_params.yaml b/flash_test/run_params.yaml new file mode 100644 index 0000000..2dc7c20 --- /dev/null +++ b/flash_test/run_params.yaml @@ -0,0 +1,56 @@ +# params to vary +init_name: Flash1D +order: 2 +mu: 1.e-3 +mesh_size: 0.0002 +current_dt: 1.e-9 +nviz: 1000 +# everything else +init_shock: False +init_flame: True +viz_level: 3 +nrestart: 1000 +mesh_angle: 0 +periodic: False +bl_ratio: 4.0 +interface_ratio: 1.0 +transfinite: True +nhealth: 100 +nstatus: 100 +t_final: 1.e-2 +mach: 4.0 +pressure_background: 101325 +temperature_background: 300 +use_species_limiter: 1 +vel_sigma: 0 +temp_sigma: 0 +nspecies: 7 +use_av: False +use_species_limiter: True +noslip: True +adiabatic: False +#mu: 1.e-5 +#kappa: 2.0 +alpha_sc: 0.05 +s0_sc: -5.0 +kappa_sc: 0.5 +integrator: compiled_lsrk54 +#integrator: euler +#inviscid_numerical_flux: hll +health_pres_min: 0.001 +health_pres_max: 1000000 +# +# wall settings +# +wall_penalty_amount: 0 +wall_time_scale: 1 +use_wall_ox: False +use_wall_mass: False +use_sponge: False +# +# filter settings +# +soln_nfilter: -1 +use_rhs_filter: False +soln_filter_cutoff: 4 +rhs_filter_cutoff: 3 diff --git a/scalability_test/data b/scalability_test/data new file mode 120000 index 0000000..62f501b --- /dev/null +++ b/scalability_test/data @@ -0,0 +1 @@ +../data/cav5_comb4/3D/scalability \ No newline at end of file diff --git a/scalability_test/driver.py b/scalability_test/driver.py new file mode 120000 index 0000000..3f406d1 --- /dev/null +++ b/scalability_test/driver.py @@ -0,0 +1 @@ +../driver.py \ No newline at end of file diff --git a/scalability_test/generate_lassen_scaling_job_script.sh b/scalability_test/generate_lassen_scaling_job_script.sh new file mode 100755 index 0000000..fa5d81f --- /dev/null +++ b/scalability_test/generate_lassen_scaling_job_script.sh @@ -0,0 +1,79 @@ +#!/bin/bash + +NONOPT_ARGS=() +while [[ $# -gt 0 ]]; do + case $1 in + -n|--nodes) + NUM_NODES="$2" + shift + shift + ;; + -e|--end) + NUM_PROCS="$2" + shift + shift + ;; + -s|--start) + NUM_PROCS_1="$2" + shift + shift + ;; + -o|--output) + OUTPUT_PATH="$2" + shift + shift + ;; + -q|--queue) + QUEUE_NAME="$2" + shift + shift + ;; + -t|--time) + TIME_LIMIT="$2" + shift + shift + ;; + -h|--help) + printf "Generates batch script for scaling tests on Lassen.\n\nUsage:\n" + printf "./generate_lassen_scaling_job_script.sh -n -s -e -o -q -t \n\n" + printf "Default: scaling_test_lassen.bsub.sh = single node scaling job to run on 1-4 GPUs on batch queue with 120m time limit.\n\n" + printf "Submit the resulting script with the \`bsub\` command.\n" + exit 1 + ;; + -*|--*) + echo "generate_lassen_scaling_job_script: Unknown option $1" + exit 1 + ;; + *) + NONOPT_ARGS+=("$1") + shift + ;; + esac +done +set -- "${NONOPT_ARGS[@]}" + + +NUM_NODES=${NUM_NODES:-"1"} +NUM_PROCS_MAX=$(( 4 * ${NUM_NODES} )) +NUM_PROCS=${NUM_PROCS:-"${NUM_PROCS_MAX}"} +NUM_PROCS_1=${NUM_PROCS_1:-"1"} +OUTPUT_PATH=${OUTPUT_PATH:-"scaling_test_lassen.bsub.sh"} +QUEUE_NAME=${QUEUE_NAME:-"pbatch"} +TIME_LIMIT=${TIME_LIMIT:-"120"} + +rm -f ${OUTPUT_PATH} +cat < ${OUTPUT_PATH} +#!/bin/bash + +#BSUB -nnodes ${NUM_NODES} +#BSUB -G uiuc +#BSUB -W ${TIME_LIMIT} +#BSUB -J scale${NUM_PROCS} +#BSUB -q ${QUEUE_NAME} +#BSUB -o scal${NUM_PROCS}.txt + +source ../emirge/config/activate_env.sh +source ../emirge/mirgecom/scripts/mirge-testing-env.sh +source ../scripts/multi_scalability.sh -p ../ -s ${NUM_PROCS_1} -n ${NUM_PROCS} + +EOF diff --git a/scalability_test/runLassenBatch.sh b/scalability_test/runLassenBatch.sh new file mode 100755 index 0000000..22dee4f --- /dev/null +++ b/scalability_test/runLassenBatch.sh @@ -0,0 +1,11 @@ +#!/bin/bash +#BSUB -nnodes 4 +#BSUB -G uiuc +#BSUB -W 180 +#BSUB -J scale16 +#BSUB -q pbatch +#BSUB -o scal16.txt + +source ../emirge/config/activate_env.sh +source ../emirge/mirgecom/scripts/mirge-testing-env.sh +source ../scripts/multi_scalability.sh -p ../ -s 1 -n 16 diff --git a/scalability_test/run_params.yaml b/scalability_test/run_params.yaml new file mode 100644 index 0000000..eca6171 --- /dev/null +++ b/scalability_test/run_params.yaml @@ -0,0 +1,69 @@ +############## +# run control +############## +nviz: 1000 +nrestart: 1000 +nhealth: -1 +nstatus: -1 +ngarbage: 10 +viz_level: 1 +viz_interval_type: 0 +#viz_interval_type: 1 +#viz_interval_type: 2 +t_viz_interval: 1.e-9 +constant_cfl: 0 +#constant_cfl: 1 +current_cfl: 0.05 +current_dt: 1.0e-13 +t_final: 2.e-11 +mesh_filename: data/actii.msh +dimen: 3 + +############## +# wall model +############## +wall_time_scale: 100 +wall_penalty_amount: 0 +wall_material: 2 + +############## +# model control +############## +use_species_limiter: True +use_av: True +alpha_sc: 0.01 +s0_sc: -5.0 +kappa_sc: 0.5 +order: 4 +# use_rhs_filter: True +# rhs_filter_cutoff: 3 + +############## +# equation of state +############## +nspecies: 7 +pyro_temp_iter: 1 +pyro_temp_tol: 1. +transport: 1 + +############## +# ignition parameters +############## +use_ignition: 2 +ignition_loc_x: 0.68 +ignition_loc_y: -0.01 +ignition_init_time: 5.e-8 +ignition_strength: 600 +ignition_duration: 1.e-8 +ignition_diameter: 0.02 +integrator: compiled_lsrk54 +#health_pres_min: 2828 +#health_pres_max: 274259 + +############## +# initialization +############## +vel_sigma: 500 +temp_sigma: 1250 +vel_sigma_inj: 5000 +temp_sigma_inj: 5000 diff --git a/scripts/multi_scalability.sh b/scripts/multi_scalability.sh new file mode 100644 index 0000000..767a1aa --- /dev/null +++ b/scripts/multi_scalability.sh @@ -0,0 +1,210 @@ +#!/bin/bash + +# Scalability driver runs a weak scaling test using a series +# of meshes of increasing sizes on [s, n] ranks. +# Example use: +# . scripts/single_scalability.sh -i -s 1 -n 16 -p ../ +# -i => optionally install the driver (default=no) +# -p => indicate path to driver top level (default=../) +# -s => number of ranks to start with (default=1) +# -n => number of ranks to end with (default=16) + +NONOPT_ARGS=() +while [[ $# -gt 0 ]]; do + case $1 in + -e|--env) + TESTING_ENV_RESOURCE="$2" + shift + shift + ;; + -c|--casename) + CASENAME_ROOT="$2" + shift + shift + ;; + -n|--nproc) + NUM_PROCS="$2" + shift + shift + ;; + -s|--start) + NUM_PROCS_1="$2" + shift + shift + ;; + -t|--testdir) + TEST_PATH="$2" + shift + shift + ;; + -o|--output) + LOG_PATH="$2" + shift + shift + ;; + -p|--path) + DRIVER_PATH="$2" + shift + shift + ;; + -i|--install) + INSTALL="YES" + shift + ;; + -*|--*) + echo "install_mirgecom: Unknown option $1" + exit 1 + ;; + *) + NONOPT_ARGS+=("$1") + shift + ;; + esac +done +set -- "${NONOPT_ARGS[@]}" + +TESTING_ENV_RESOURCE=${TESTING_ENV_RESOURCE:-""} +CASENAME_ROOT=${CASENAME_ROOT:-""} +LOG_PATH=${LOG_PATH:-"log_data"} +DRIVER_PATH=${DRIVER_PATH:-"."} +NUM_PROCS=${NUM_PROCS:-"4"} +NUM_PROCS_1=${NUM_PROCS_1:-"1"} +INSTALL=${INSTALL:-"NO"} +TEST_PATH=${TEST_PATH:-"scalability_test"} + +if [[ ! -z "${CASENAME_ROOT}" ]]; then + printf "Casename file prefix: ${CASENAME_ROOT}\n" +fi + +if [[ ! -z "${TESTING_ENV_RESOURCE}" ]]; then + printf "Testing environment resource file: ${TESTING_ENV_RESOURCE}\n" + source ${TESTING_ENV_RESOURCE} +fi + +mkdir -p ${LOG_PATH} +cd ${LOG_PATH} +LOG_PATH=$(pwd) +cd - + +# Set defaults for these in case they didn't get +# set by the resource file. +# At this level, MPI_EXEC, and PARALLEL_SPAWNER are all that +# are needed to make this script platform-agnostic +MPI_EXEC=${MIRGE_MPI_EXEC:-"mpiexec"} +PARALLEL_SPAWNER=${MIRGE_PARALLEL_SPAWNER:-""} + +# path management junk +# Set path to top of driver (assume cwd if unset) +origin=$(pwd) + +# Driver scripts need to be run from their home directory to find +# their own resources. +cd ${DRIVER_PATH} +DRIVER_PATH=$(pwd) + +printf "Driver directory: ${DRIVER_PATH}\n" + +if [ "${INSTALL}" == "YES" ]; then + printf "Installing driver...\n" + + python -m pip install -e . + + printf "Driver installed.\n" +fi + +date + +# Tracks success/fail of multiple runs +# (i.e. one for each mesh size) +declare -i numfail=0 +declare -i numsuccess=0 +succeeded_tests="" +failed_tests="" +if [[ ! -z ${CASENAME_ROOT} ]];then + CASENAME_ROOT="${CASENAME_ROOT}_" +fi + +# - The casename thing is a "nice to have" as it allows some control +# to caller for what the sqlite files are named. +printf "Running parallel timing tests...\n" +test_path=${TEST_PATH} +test_name="prediction-scalability" +printf "* Running ${test_name} test in ${test_path}.\n" +running_casename_base="${CASENAME_ROOT}${test_name}" + +cd ${DRIVER_PATH}/${test_path} + +nrank=${NUM_PROCS_1} +while [ $nrank -le $NUM_PROCS ]; do + + if [[ "${nrank}" == "1" ]]; then + msize="48" + nelem="24036" + elif [[ "${nrank}" == "2" ]]; then + msize="30.5" + nelem="47908" + elif [[ "${nrank}" == "4" ]]; then + msize="21.5" + nelem="96425" + elif [[ "${nrank}" == "8" ]]; then + msize="16.05" + nelem="192658" + elif [[ "${nrank}" == "16" ]]; then + msize="12.3" + nelem="383014" + elif [[ "${nrank}" == "32" ]]; then + msize="9.535" + nelem="768891" + elif [[ "${nrank}" == "64" ]]; then + msize="7.45" + nelem="1536571" + elif [[ "${nrank}" == "128" ]]; then + msize="5.8545" + nelem="3071835" + elif [[ "${nrank}" == "256" ]]; then + msize="4.61" + nelem="6142150" + elif [[ "${nrank}" == "512" ]]; then + msize="3.6389" + nelem="12283032" + fi + + + casename="${running_casename_base}_np${nrank}" + printf "** Running ${casename} on ${nrank} ranks.\n" + + cd data + rm actii.msh + ./mkmsh --size=${msize} --nelem=${nelem} --link + cd ../ + + set -x + $MPI_EXEC -n ${nrank} $PARALLEL_SPAWNER python -u -m mpi4py driver.py -c ${casename} -g ${LOG_PATH} -i run_params.yaml --log --lazy + return_code=$? + set +x + + mv viz_data viz_data_${nrank} + mv restart_data restart_data_${nrank} + + if [[ $return_code -eq 0 ]] + then + ((numsuccess=numsuccess+1)) + echo "** ${test_path}/${casename} $succeeded." + succeeded_tests="$succeeded_tests ${test_path}/${casename}" + else + ((numfail=numfail+1)) + echo "** ${test_path}/${casename} failed." + failed_tests="$failed_tests ${test_path}/${casename}" + fi + + nrank=$(( $nrank * 2 )) + +done + +date +printf "Scaling/timing tests done.\n" +printf "Passing tests: ${succeeded_tests}\n" +printf "Failing tests: ${failed_tests}\n" + +cd ${origin} +return $numfail diff --git a/scripts/single_scalability.sh b/scripts/single_scalability.sh new file mode 100644 index 0000000..92edda2 --- /dev/null +++ b/scripts/single_scalability.sh @@ -0,0 +1,178 @@ +#!/bin/bash + +# Runs the single node scalability (1, 2, 4) ranks +# Example use: +# . scripts/single_scalability.sh -i -p ../ +# -i => optionally install the driver +# -p => indicate the path to top level of driver + +NONOPT_ARGS=() +while [[ $# -gt 0 ]]; do + case $1 in + -e|--env) + TESTING_ENV_RESOURCE="$2" + shift + shift + ;; + -c|--casename) + CASENAME_ROOT="$2" + shift + shift + ;; + -i|--install) + INSTALL="YES" + shift + ;; + -o|--output) + LOG_PATH="$2" + shift + shift + ;; + -p|--path) + DRIVER_PATH="$2" + shift + shift + ;; + -t|--testdir) + TEST_PATH="$2" + shift + shift + ;; + -*|--*) + echo "install_mirgecom: Unknown option $1" + exit 1 + ;; + *) + NONOPT_ARGS+=("$1") + shift + ;; + esac +done +set -- "${NONOPT_ARGS[@]}" + +TESTING_ENV_RESOURCE=${TESTING_ENV_RESOURCE:-""} +CASENAME_ROOT=${CASENAME_ROOT:-""} +LOG_PATH=${LOG_PATH:-"log_data"} +DRIVER_PATH=${DRIVER_PATH:-"."} +INSTALL=${INSTALL:-"NO"} +TEST_PATH=${TEST_PATH:-"scalability_test"} + +if [[ ! -z "${CASENAME_ROOT}" ]]; then + printf "Casename file prefix: ${CASENAME_ROOT}\n" +fi + +if [[ ! -z "${TESTING_ENV_RESOURCE}" ]]; then + printf "Testing environment resource file: ${TESTING_ENV_RESOURCE}\n" + source ${TESTING_ENV_RESOURCE} +fi + +mkdir -p ${LOG_PATH} +cd ${LOG_PATH} +LOG_PATH=$(pwd) +cd - + +# Set defaults for these in case they didn't get +# set by the resource file. +# At this level, MPI_EXEC, and PARALLEL_SPAWNER are all that +# are needed to make this script platform-agnostic +MPI_EXEC=${MIRGE_MPI_EXEC:-"mpiexec"} +PARALLEL_SPAWNER=${MIRGE_PARALLEL_SPAWNER:-""} + +# path management junk +# Set path to top of driver (assume cwd if unset) +origin=$(pwd) + +# Driver scripts need to be run from their home directory to find +# their own resources. +cd ${DRIVER_PATH} +DRIVER_PATH=$(pwd) + +printf "Driver directory: ${DRIVER_PATH}\n" + +if [ "${INSTALL}" == "YES" ]; then + printf "Installing driver...\n" + + python -m pip install -e . + + printf "Driver installed.\n" +fi + +date + +declare -i numfail=0 +declare -i numsuccess=0 +succeeded_tests="" +failed_tests="" +if [[ ! -z ${CASENAME_ROOT} ]];then + CASENAME_ROOT="${CASENAME_ROOT}_" +fi + +# - The casename thing is a "nice to have" as it allows some control +# to caller for what the sqlite files are named. +# +printf "Running parallel timing tests...\n" + +test_path=${TEST_PATH} +test_name="prediction-scalability" + +printf "* Running ${test_name} test in ${test_path}.\n" +running_casename_base="${CASENAME_ROOT}${test_name}" + +cd ${DRIVER_PATH}/${test_path} + +for nrank in 1 2 4; do + + if [[ "${nrank}" == "1" ]]; then + msize="48" + nelem="24036" + elif [[ "${nrank}" == "2" ]]; then + msize="30.5" + nelem="47908" + elif [[ "${nrank}" == "4" ]]; then + msize="21.5" + nelem="96425" + elif [[ "${nrank}" == "8" ]]; then + msize="16.05" + nelem="192658" + elif [[ "${nrank}" == "16" ]]; then + msize="12.3" + nelem="383014" + fi + + casename="${running_casename_base}_np${nrank}" + printf "** Running ${casename} on ${nrank} ranks.\n" + + cd data + rm actii.msh + ./mkmsh --size=${msize} --nelem=${nelem} --link + cd ../ + + set -x + $MPI_EXEC -n ${nrank} $PARALLEL_SPAWNER python -u -m mpi4py driver.py -c ${casename} -g ${LOG_PATH} -i run_params.yaml --log --lazy + return_code=$? + set +x + + mv viz_data viz_data_${nrank} + mv restart_data restart_data_${nrank} + + + if [[ $return_code -eq 0 ]] + then + ((numsuccess=numsuccess+1)) + echo "** ${test_path}/${casename} $succeeded." + succeeded_tests="$succeeded_tests ${test_path}/${casename}" + else + ((numfail=numfail+1)) + echo "** ${test_path}/${casename} failed." + failed_tests="$failed_tests ${test_path}/${casename}" + fi + +done + +date +printf "Scaling/timing tests done.\n" +printf "Passing tests: ${succeeded_tests}\n" +printf "Failing tests: ${failed_tests}\n" + +cd ${origin} +return $numfail diff --git a/scripts/smoke_test.sh b/scripts/smoke_test.sh index 6b9214e..3ce4859 100755 --- a/scripts/smoke_test.sh +++ b/scripts/smoke_test.sh @@ -38,7 +38,7 @@ cd ${TOP_PATH} TOP_PATH=$(pwd) printf "Driver directory: ${TOP_PATH}\n" -python -m pip install . +python -m pip install -e . date @@ -60,12 +60,7 @@ do # Create 3d mesh if not already there if [[ "${test_name}" == *"_3d"* ]]; then cd data - rm actii.msh - if [[ -f "actii_24110.msh" ]]; then - ln -s actii_24110.msh actii.msh - else - ./mkmsh --size=48 --link # will not overwrite if exists - fi + ./mkmsh --size=48 --nelem=24110 --link # will not overwrite if exists cd ../ fi @@ -99,6 +94,11 @@ do cd ${TOP_PATH}/${test_path} # Create 3d mesh unless already there +<<<<<<< HEAD + cd data + ./mkmsh --size=30.5 --nelem=47908 --link # will not overwrite if it exists + cd ../ +======= if [ "${test_name}" == "smoke_test_ks_3d" ]; then cd data rm actii.msh @@ -109,6 +109,7 @@ do fi cd ../ fi +>>>>>>> main # Run the case $MPI_EXEC -n 2 $PARALLEL_SPAWNER python -u -m mpi4py driver.py -i run_params.yaml --log --lazy diff --git a/y3prediction/actii_y3.py b/y3prediction/actii_y3.py index 4a1411e..f755efe 100644 --- a/y3prediction/actii_y3.py +++ b/y3prediction/actii_y3.py @@ -667,3 +667,141 @@ def __call__(self, dcoll, x_vec, eos, *, time=0.0): energy=energy, species_mass=mass*y ) + + +def get_mesh_data(mesh_filename=None, dim=None): + if any([mesh_filename, dim]) is None: + raise ValueError("ACTII grid initialization requires filename and" + " dimension.") + + from meshmode.mesh.io import read_gmsh + print(f"Reading {dim}D mesh from {mesh_filename}") + + mesh, tag_to_elements = read_gmsh( + mesh_filename, force_ambient_dim=dim, + return_tag_to_elements_map=True) + + volume_to_tags = { + "fluid": ["fluid"], + "wall": ["wall_insert", "wall_surround"]} + + return mesh, tag_to_elements, volume_to_tags + + +def get_boundaries(dcoll, actx, dd_vol_fluid, dd_vol_wall, use_injection, + quadrature_tag, gas_model, noslip, adiabatic, + temp_wall, target_fluid_state): + + from mirgecom.boundary import ( + PrescribedFluidBoundary, + IsothermalWallBoundary, + AdiabaticSlipBoundary, + AdiabaticNoslipWallBoundary, + DummyBoundary + ) + from mirgecom.diffusion import ( + DirichletDiffusionBoundary + ) + from mirgecom.simutil import force_evaluation + from mirgecom.gas_model import project_fluid_state + + # use dummy boundaries to setup the smoothness state for the target + wall_bnd = dd_vol_fluid.trace("isothermal_wall") + inflow_bnd = dd_vol_fluid.trace("inflow") + outflow_bnd = dd_vol_fluid.trace("outflow") + inj_bnd = dd_vol_fluid.trace("injection") + flow_bnd = dd_vol_fluid.trace("flow") + wall_ffld_bnd = dd_vol_wall.trace("wall_farfield") + + if use_injection: + target_boundaries = { + flow_bnd.domain_tag: # pylint: disable=no-member + DummyBoundary(), + wall_bnd.domain_tag: # pylint: disable=no-member + IsothermalWallBoundary() + } + else: + target_boundaries = { + inflow_bnd.domain_tag: # pylint: disable=no-member + DummyBoundary(), + outflow_bnd.domain_tag: # pylint: disable=no-member + DummyBoundary(), + inj_bnd.domain_tag: # pylint: disable=no-member + IsothermalWallBoundary(), + wall_bnd.domain_tag: # pylint: disable=no-member + IsothermalWallBoundary() + } + + ################################## + # Set up the boundary conditions # + ################################## + + def get_target_state_on_boundary(btag): + return project_fluid_state( + dcoll, dd_vol_fluid, + dd_vol_fluid.trace(btag).with_discr_tag(quadrature_tag), + target_fluid_state, gas_model + ) + + flow_ref_state = \ + get_target_state_on_boundary("flow") + + flow_ref_state = force_evaluation(actx, flow_ref_state) + + def _target_flow_state_func(**kwargs): + return flow_ref_state + + flow_boundary = PrescribedFluidBoundary( + boundary_state_func=_target_flow_state_func) + + inflow_ref_state = \ + get_target_state_on_boundary("inflow") + + inflow_ref_state = force_evaluation(actx, inflow_ref_state) + + def _target_inflow_state_func(**kwargs): + return inflow_ref_state + + inflow_boundary = PrescribedFluidBoundary( + boundary_state_func=_target_inflow_state_func) + + outflow_ref_state = \ + get_target_state_on_boundary("outflow") + + outflow_ref_state = force_evaluation(actx, outflow_ref_state) + + def _target_outflow_state_func(**kwargs): + return outflow_ref_state + + outflow_boundary = PrescribedFluidBoundary( + boundary_state_func=_target_outflow_state_func) + #outflow_pressure = 2000 + #outflow_boundary = PressureOutflowBoundary(outflow_pressure) + + if noslip: + if adiabatic: + fluid_wall = AdiabaticNoslipWallBoundary() + else: + fluid_wall = IsothermalWallBoundary(temp_wall) + else: + fluid_wall = AdiabaticSlipBoundary() + + wall_farfield = DirichletDiffusionBoundary(temp_wall) + + if use_injection: + fluid_boundaries = { + flow_bnd.domain_tag: flow_boundary, # pylint: disable=no-member + wall_bnd.domain_tag: fluid_wall # pylint: disable=no-member + } + else: + fluid_boundaries = { + inflow_bnd.domain_tag: inflow_boundary, # pylint: disable=no-member + outflow_bnd.domain_tag: outflow_boundary, # pylint: disable=no-member + inj_bnd.domain_tag: fluid_wall, # pylint: disable=no-member + wall_bnd.domain_tag: fluid_wall # pylint: disable=no-member + } + + wall_boundaries = { + wall_ffld_bnd.domain_tag: wall_farfield # pylint: disable=no-member + } + return fluid_boundaries, wall_boundaries, target_boundaries diff --git a/y3prediction/flash_utils.py b/y3prediction/flash_utils.py new file mode 100644 index 0000000..8d7490f --- /dev/null +++ b/y3prediction/flash_utils.py @@ -0,0 +1,1027 @@ +import numpy as np + + +def setup_flame(mechanism_name="uiuc", fuel_name="C2H4", + temperature_unburned=300.0, + stoich_ratio=None, equiv_ratio=1.0, ox_di_ratio=.21, + oxidizer_name="O2", inert_name="N2", + pressure_unburned=None): + # {{{ Set up burned/unburned state using Cantera + + import cantera + from mirgecom.mechanisms import get_mechanism_input + mech_input = get_mechanism_input(mechanism_name) + cantera_soln = cantera.Solution(name="gas", yaml=mech_input) + nspecies = cantera_soln.n_species + species_names = cantera_soln.species_names + + # Initial temperature, pressure, and mixutre mole fractions are needed to + # set up the initial state in Cantera. + # Parameters for calculating the amounts of fuel, oxidizer, and inert species + if stoich_ratio is None: + if fuel_name == "C2H4": + stoich_ratio = 3.0 + elif fuel_name == "H2": + stoich_ratio = 0.5 + else: + stoich_ratio = 1.0 + + # Grab the array indices for the specific species, ethylene, oxygen, and nitrogen + i_fu = cantera_soln.species_index(fuel_name) + i_ox = cantera_soln.species_index(oxidizer_name) + i_di = cantera_soln.species_index(inert_name) + x = np.zeros(nspecies) + + # Set the species mole fractions according to our desired fuel/air mixture + x[i_fu] = (ox_di_ratio*equiv_ratio)/(stoich_ratio+ox_di_ratio*equiv_ratio) + x[i_ox] = stoich_ratio*x[i_fu]/equiv_ratio + x[i_di] = (1.0-ox_di_ratio)*x[i_ox]/ox_di_ratio + + one_atm = cantera.one_atm # pylint: disable=no-member + pres_unburned = one_atm if pressure_unburned is None else pressure_unburned + + # Let the user know about how Cantera is being initilized + print(f"Input state (T,P,X) = ({temperature_unburned}, {pres_unburned}, {x}") + # Set Cantera internal gas temperature, pressure, and mole fractios + cantera_soln.TPX = temperature_unburned, pres_unburned, x + # Pull temperature, total density, mass fractions, and pressure from Cantera + # We need total density, and mass fractions to initialize the fluid/gas state. + y_unburned = np.zeros(nspecies) + can_t, rho_unburned, y_unburned = cantera_soln.TDY + + # now find the conditions for the burned gas + cantera_soln.equilibrate("TP") + temp_burned, rho_burned, y_burned = cantera_soln.TDY + pres_burned = cantera_soln.P + + print(f"Mechanism species names {species_names}") + print(f"Unburned (T,P,Y) = ({temperature_unburned}, {pres_unburned}, " + f"{y_unburned}") + print(f"Burned (T,P,Y) = ({temp_burned}, {pres_burned}, {y_burned}") + + return {"unburned": (pres_unburned, temperature_unburned, y_unburned, + rho_unburned), + "burned": (pres_burned, temp_burned, y_burned, rho_burned)} + + +class Flash1D: + r"""Solution initializer for flow with a shock and/or flame. + + This initializer creates a physics-consistent flow solution + given an initial thermal state (pressure, temperature) for the pre and post + burned and shocked regions and an EOS. + + The solution varies across a planar shock/flame interfaces defined by a tanh + function located at (shock,flame)_location for pressure, temperature, velocity, + and mass fraction + + .. automethod:: __init__ + .. automethod:: __call__ + """ + + def __init__( + self, *, dim=3, nspecies=0, + shock_normal_dir=None, flame_normal_dir=None, + shock_location=None, flame_location=None, + temperature_unshocked=None, temperature_shocked=None, + temperature_unburned=None, temperature_burned=None, + pressure_unshocked=None, pressure_shocked=None, + pressure_unburned=None, pressure_burned=None, + species_mass_fractions_shocked=None, + species_mass_fractions_unshocked=None, + species_mass_fractions_burned=None, + species_mass_fractions_unburned=None, + velocity_unshocked=None, velocity_shocked=None, + velocity_cross=None, + convective_velocity=None, sigma=0.5, + temp_sigma=0., vel_sigma=0., temp_wall=300. + ): + r"""Initialize mixture parameters. + + Parameters + ---------- + dim: int + specifies the number of dimensions for the solution + nspecies: int + specifies the number of mixture species + shock_normal_dir: numpy.ndarray + specifies the orientation of the discontinuity + shock_location: numpy.ndarray + init/fixed location of discontinuity + flame_normal_dir: numpy.ndarray + specifies the orientation of the flame front + flame_location: numpy.ndarray + fixed location of the flame front + pressure_shocked: float + pressure in the shocked region + pressure_unshocked: float + pressure in the unshocked region + temperature_shocked: float + temperature in the shocked region + temperature_unshocked: float + temperature in the unshocked region + velocity_shocked: numpy.ndarray + velocity (vector) in the shocked region + velocity_unshocked: numpy.ndarray + velocity (vector) in the unshocked region + species_mass_fractions_shocked: numpy.ndarray + species mass fractions in the shocked region + species_mass_fractions_unshocked: numpy.ndarray + species mass fractions in the unshocked region + pressure_burned: float + pressure in the shocked region + pressure_unburned: float + pressure in the unshocked region + temperature_burned: float + temperature in the shocked region + temperature_unburned: float + temperature in the unshocked region + species_mass_fractions_burned: numpy.ndarray + species mass fractions in the shocked region + species_mass_fractions_unburned: numpy.ndarray + species mass fractions in the unshocked region + sigma: float + sharpness parameter for shock and flame fronts + velocity_cross: numpy.ndarray + velocity (vector) tangent to the shock + temp_sigma: float + near-wall temperature relaxation parameter + vel_sigma: float + near-wall velocity relaxation parameter + """ + self.do_shock = shock_location is not None + self.do_flame = flame_location is not None + if not self.do_flame and not self.do_shock: + raise ValueError("Invalid arguments for flame and shock init.") + + if self.do_flame: + if nspecies == 0: + raise ValueError("nspecies must be set for flame initialization.") + if any([pressure_burned, pressure_unburned, temperature_burned, + temperature_unburned, species_mass_fractions_burned, + species_mass_fractions_unburned, flame_normal_dir]) is None: + raise ValueError("burned/unburned states underspecified for" + " flame init.") + + if self.do_shock: + if any([pressure_shocked, pressure_unshocked, temperature_shocked, + temperature_unshocked, shock_normal_dir]) is None: + raise ValueError("shocked/unshocked states underspecified for " + "flame init.") + if not self.do_flame and nspecies > 0: + if any([species_mass_fractions_shocked, + species_mass_fractions_unshocked]) is None: + raise ValueError("y shocked/unshocked must be given for" + " multispecies shock.") + + if velocity_shocked is None: + velocity_shocked = np.zeros(shape=(dim,)) + if velocity_unshocked is None: + velocity_unshocked = np.zeros(shape=(dim,)) + if velocity_cross is None: + velocity_cross = np.zeros(shape=(dim,)) + if convective_velocity is None: + convective_velocity = np.zeros(shape=(dim,)) + + self._nspecies = nspecies + self._dim = dim + self._sigma = sigma + self._temp_sigma = temp_sigma + self._vel_sigma = vel_sigma + self._temp_wall = temp_wall + + self._shock_loc = shock_location + self._flame_loc = flame_location + self._shock_normal = shock_normal_dir + self._flame_normal = flame_normal_dir + + self._u_shocked = velocity_shocked + self._u_unshocked = velocity_unshocked + self._ut = velocity_cross + self._uc = convective_velocity + self._p_shocked = pressure_shocked + self._p_unshocked = pressure_unshocked + self._t_shocked = temperature_shocked + self._t_unshocked = temperature_unshocked + + self._t_burned = temperature_burned + self._t_unburned = temperature_unburned + self._p_burned = pressure_burned + self._p_unburned = pressure_unburned + self._y_burned = species_mass_fractions_burned + self._y_unburned = species_mass_fractions_unburned + self._y_shocked = species_mass_fractions_shocked + self._y_unshocked = species_mass_fractions_unshocked + + def __call__(self, x_vec, eos, *, time=0.0): + """Create the mixture state at locations *x_vec*. + + Parameters + ---------- + x_vec: numpy.ndarray + Coordinates at which solution is desired + eos: + Mixture-compatible equation-of-state object must provide + these functions: + `eos.get_density` + `eos.get_internal_energy` + time: float + Time at which solution is desired. The location is (optionally) + dependent on time + """ + if x_vec.shape != (self._dim,): + raise ValueError(f"Position vector has unexpected dimensionality," + f" expected {self._dim}.") + + xpos = x_vec[0] + ypos = x_vec[1] + ones = 0*xpos + 1. + # if self._dim == 3: + # zpos = x_vec[2] + + actx = xpos.array_context + + # Set up flame states + flame_rd = 1.*ones # unburned everywhere + if self.do_flame: + flame_x0 = self._flame_loc + flame_rd = np.dot(x_vec - flame_x0, self._flame_normal) / self._sigma + dp = self._p_unburned - self._p_burned + dy = self._y_unburned - self._y_burned + dt = self._t_unburned - self._t_burned + weight = 0.5*(1.0 + actx.np.tanh(flame_rd)) + y_flame = self._y_burned + dy*weight + p_flame = self._p_burned + dp*weight + t_flame = self._t_burned + dt*weight + pressure = p_flame + temperature = t_flame + y = y_flame + + # Set up shock states + if self.do_shock: + shock_x0 = self._shock_loc + shock_rd = np.dot(x_vec - shock_x0, self._shock_normal) / self._sigma + dp = self._p_unshocked - self._p_shocked + dt = self._t_unshocked - self._t_shocked + weight = 0.5*(1.0 + actx.np.tanh(shock_rd)) + p_shock = self._p_shocked + dp*weight + t_shock = self._t_shocked + dt*weight + pressure = p_shock + temperature = t_shock + if self._nspecies > 0: + dy = self._y_unshocked - self._y_shocked + y = self._y_shocked + dy*weight + else: + y = None + + velocity = \ + self._u_shocked + (self._u_unshocked-self._u_shocked)*weight + self._ut + + if self.do_shock and self.do_flame: + # Use Y from flame + # Use pressure from: + # - shock (if shocked region) + # - flame (if unshocked region) + # Use max temperature + y = y_flame + pressure = actx.np.where(actx.np.less(shock_rd, 0), p_shock, p_flame) + temperature = actx.np.max(t_shock, t_flame) + + # modify the temperature in the near wall region to match the + # isothermal boundaries + y_top = 0.01 + y_bottom = -0.01 + if self._temp_sigma > 0: + sigma = self._temp_sigma + wall_temperature = self._temp_wall + smoothing_top = smooth_step(actx, -sigma*(ypos - y_top)) + smoothing_bottom = smooth_step(actx, sigma*(ypos - y_bottom)) + temperature = \ + (wall_temperature + + (temperature - wall_temperature)*smoothing_top*smoothing_bottom) + + # modify the velocity in the near wall region to match the + # noslip boundaries + sigma = self._vel_sigma + if sigma > 0: + smoothing_top = smooth_step(actx, -sigma*(ypos - y_top)) + smoothing_bottom = smooth_step(actx, sigma*(ypos - y_bottom)) + velocity[0] = velocity[0]*smoothing_top*smoothing_bottom + + mass = eos.get_density(pressure, temperature, species_mass_fractions=y) + specmass = mass * y if y is not None else None + mom = mass * velocity + internal_energy = eos.get_internal_energy(temperature, + species_mass_fractions=y) + + kinetic_energy = 0.5 * np.dot(velocity, velocity) + energy = mass * (internal_energy + kinetic_energy) + + from mirgecom.fluid import make_conserved + return make_conserved(dim=self._dim, mass=mass, energy=energy, + momentum=mom, species_mass=specmass) + + +def get_flash_mesh_2d(size, bl_ratio, interface_ratio, angle=0., + transfinite=False, + left_boundary_loc=0., + bottom_boundary_loc=-0.01, + aft_boundary_loc=-0.01): + """Generate a grid using `gmsh`. + + """ + + dim = 2 + height = 0.02 + fluid_length = 0.1 + wall_length = 0.02 + bottom_inflow = np.zeros(shape=(dim,)) + top_inflow = np.zeros(shape=(dim,)) + bottom_interface = np.zeros(shape=(dim,)) + top_interface = np.zeros(shape=(dim,)) + bottom_wall = np.zeros(shape=(dim,)) + top_wall = np.zeros(shape=(dim,)) + + # rotate the mesh around the bottom-left corner + theta = angle/180.*np.pi/2. + bottom_inflow = [left_boundary_loc, bottom_boundary_loc, aft_boundary_loc] + top_inflow[0] = bottom_inflow[0] - height*np.sin(theta) + top_inflow[1] = bottom_inflow[1] + height*np.cos(theta) + + bottom_interface[0] = bottom_inflow[0] + fluid_length*np.cos(theta) + bottom_interface[1] = bottom_inflow[1] + fluid_length*np.sin(theta) + top_interface[0] = top_inflow[0] + fluid_length*np.cos(theta) + top_interface[1] = top_inflow[1] + fluid_length*np.sin(theta) + + bottom_wall[0] = bottom_interface[0] + wall_length*np.cos(theta) + bottom_wall[1] = bottom_interface[1] + wall_length*np.sin(theta) + top_wall[0] = top_interface[0] + wall_length*np.cos(theta) + top_wall[1] = top_interface[1] + wall_length*np.sin(theta) + + from meshmode.mesh.io import ( + generate_gmsh, + ScriptSource + ) + + # for 2D, the line segments/surfaces need to be specified clockwise to + # get the correct facing (right-handed) surface normals + # for 3D this is the back plane + my_string = \ + f""" + Point(1) = {{ {bottom_inflow[0]}, {bottom_inflow[1]}, + {bottom_inflow[2]}, {size}}}; + Point(2) = {{ {bottom_interface[0]}, {bottom_interface[1]}, + {bottom_interface[2]}, {size}}}; + Point(3) = {{ {top_interface[0]}, {top_interface[1]}, + {top_interface[2]}, {size}}}; + Point(4) = {{ {top_inflow[0]}, {top_inflow[1]}, + {top_inflow[2]}, {size}}}; + Point(5) = {{ {bottom_wall[0]}, {bottom_wall[1]}, + {bottom_wall[2]}, {size}}}; + Point(6) = {{ {top_wall[0]}, {top_wall[1]}, + {top_wall[2]}, {size}}}; + Line(1) = {{1, 2}}; + Line(2) = {{2, 3}}; + Line(3) = {{3, 4}}; + Line(4) = {{4, 1}}; + Line(5) = {{3, 6}}; + Line(6) = {{2, 5}}; + Line(7) = {{5, 6}}; + Line Loop(1) = {{-4, -3, -2, -1}}; + Line Loop(2) = {{2, 5, -7, -6}}; + Plane Surface(1) = {{1}}; + Plane Surface(2) = {{2}}; + + Physical Surface('fluid') = {{1}}; + Physical Surface('wall') = {{2}}; + Physical Curve('fluid_inflow') = {{4}}; + Physical Curve('fluid_wall') = {{1,3}}; + Physical Curve('fluid_wall_top') = {{3}}; + Physical Curve('fluid_wall_bottom') = {{1}}; + Physical Curve('interface') = {{2}}; + Physical Curve('wall_farfield') = {{5, 6, 7}}; + Physical Curve('solid_wall_top') = {{5}}; + Physical Curve('solid_wall_bottom') = {{6}}; + Physical Curve('solid_wall_end') = {{7}}; + """ + + if transfinite: + transfinite_string = \ + f""" + Transfinite Curve {{1, 3}} = {0.1} / {size}; + Transfinite Curve {{5, 6}} = {0.02} / {size}; + Transfinite Curve {{-2, 4, 7}}={0.02}/{size} Using Bump 1/{bl_ratio}; + Transfinite Surface {{1, 2}} Right; + Mesh.MeshSizeExtendFromBoundary = 0; + Mesh.MeshSizeFromPoints = 0; + Mesh.MeshSizeFromCurvature = 0; + Mesh.Algorithm = 5; + Mesh.OptimizeNetgen = 1; + Mesh.Smoothing = 0; + """ + my_string = my_string + transfinite_string + else: + mesh_tail = \ + f""" + // Create distance field from curves, excludes cavity + Field[1] = Distance; + Field[1].CurvesList = {{1,3}}; + Field[1].NumPointsPerCurve = 100000; + + //Create threshold field that varrries element size near boundaries + Field[2] = Threshold; + Field[2].InField = 1; + Field[2].SizeMin = {size} / {bl_ratio}; + Field[2].SizeMax = {size}; + Field[2].DistMin = 0.0002; + Field[2].DistMax = 0.005; + Field[2].StopAtDistMax = 1; + + // background mesh size + Field[3] = Box; + Field[3].XMin = 0.; + Field[3].XMax = 1.0; + Field[3].YMin = -1.0; + Field[3].YMax = 1.0; + Field[3].VIn = {size}; + + // Create distance field from curves, excludes cavity + Field[4] = Distance; + Field[4].CurvesList = {{2}}; + Field[4].NumPointsPerCurve = 100000; + + //Create threshold field that varrries element size near boundaries + Field[5] = Threshold; + Field[5].InField = 4; + Field[5].SizeMin = {size} / {interface_ratio}; + Field[5].SizeMax = {size}; + Field[5].DistMin = 0.0002; + Field[5].DistMax = 0.005; + Field[5].StopAtDistMax = 1; + + // take the minimum of all defined meshing fields + Field[100] = Min; + Field[100].FieldsList = {{2, 3, 5}}; + Background Field = 100; + + Mesh.MeshSizeExtendFromBoundary = 0; + Mesh.MeshSizeFromPoints = 0; + Mesh.MeshSizeFromCurvature = 0; + + Mesh.Algorithm = 5; + Mesh.OptimizeNetgen = 1; + Mesh.Smoothing = 100; + """ + my_string = my_string + mesh_tail + + #print(my_string) + from functools import partial + generate_mesh = partial(generate_gmsh, ScriptSource(my_string, "geo"), + force_ambient_dim=dim, dimensions=dim, + target_unit="M", return_tag_to_elements_map=True) + + return generate_mesh + + +def get_flash_mesh_3d(size, bl_ratio, interface_ratio, + angle=0., transfinite=False, + left_boundary_loc=-1, right_boundary_loc=1, + bottom_boundary_loc=-1, top_boundary_loc=1): + """Generate a grid using `gmsh`. + + """ + + dim = 3 + height = 0.02 + fluid_length = 0.1 + wall_length = 0.02 + bottom_inflow = np.zeros(shape=(dim,)) + top_inflow = np.zeros(shape=(dim,)) + bottom_interface = np.zeros(shape=(dim,)) + top_interface = np.zeros(shape=(dim,)) + bottom_wall = np.zeros(shape=(dim,)) + top_wall = np.zeros(shape=(dim,)) + + # rotate the mesh around the bottom-left corner + theta = angle/180.*np.pi/2. + bottom_inflow[0] = 0.0 + bottom_inflow[1] = -0.01 + top_inflow[0] = bottom_inflow[0] - height*np.sin(theta) + top_inflow[1] = bottom_inflow[1] + height*np.cos(theta) + + bottom_interface[0] = bottom_inflow[0] + fluid_length*np.cos(theta) + bottom_interface[1] = bottom_inflow[1] + fluid_length*np.sin(theta) + top_interface[0] = top_inflow[0] + fluid_length*np.cos(theta) + top_interface[1] = top_inflow[1] + fluid_length*np.sin(theta) + + bottom_wall[0] = bottom_interface[0] + wall_length*np.cos(theta) + bottom_wall[1] = bottom_interface[1] + wall_length*np.sin(theta) + top_wall[0] = top_interface[0] + wall_length*np.cos(theta) + top_wall[1] = top_interface[1] + wall_length*np.sin(theta) + + from meshmode.mesh.io import ( + generate_gmsh, + ScriptSource + ) + + # for 2D, the line segments/surfaces need to be specified clockwise to + # get the correct facing (right-handed) surface normals + my_string = \ + f""" + Point(1) = {{ {bottom_inflow[0]}, {bottom_inflow[1]}, 0, {size}}}; + Point(2) = {{ {bottom_interface[0]}, {bottom_interface[1]}, 0, {size}}}; + Point(3) = {{ {top_interface[0]}, {top_interface[1]}, 0, {size}}}; + Point(4) = {{ {top_inflow[0]}, {top_inflow[1]}, 0, {size}}}; + Point(5) = {{ {bottom_wall[0]}, {bottom_wall[1]}, 0, {size}}}; + Point(6) = {{ {top_wall[0]}, {top_wall[1]}, 0, {size}}}; + Line(1) = {{1, 2}}; + Line(2) = {{2, 3}}; + Line(3) = {{3, 4}}; + Line(4) = {{4, 1}}; + Line(5) = {{3, 6}}; + Line(6) = {{2, 5}}; + Line(7) = {{5, 6}}; + Line Loop(1) = {{-4, -3, -2, -1}}; + Line Loop(2) = {{2, 5, -7, -6}}; + Plane Surface(1) = {{1}}; + Plane Surface(2) = {{2}}; + + fluid_surf_vec[] = Extrude {{0, 0, {height} }} {{ Surface{{1}}; }}; + wall_surf_vec[] = Extrude {{0, 0, {height} }} {{ Surface{{2}}; }}; + + Coherence; + + Physical Volume('fluid') = {{1}}; + Physical Volume('wall') = {{2}}; + Physical Surface('fluid_inflow') = {{16}}; + Physical Surface('fluid_wall') = {{1, 20, 28, 29}}; + Physical Surface('fluid_wall_top') = {{20}}; + Physical Surface('fluid_wall_bottom') = {{28}}; + Physical Surface('fluid_wall_fore') = {{29}}; + Physical Surface('fluid_wall_aft') = {{1}}; + Physical Surface('wall_farfield') = {{2, 42, 46, 50, 51}}; + Physical Surface('solid_wall_top') = {{42}}; + Physical Surface('solid_wall_bottom') = {{50}}; + Physical Surface('solid_wall_fore') = {{51}}; + Physical Surface('solid_wall_aft') = {{2}}; + Physical Surface('solid_wall_end') = {{46}}; + """ + + if transfinite: + # MJA, fixme + transfinite_string = \ + f""" + Transfinite Curve {{1, 3}} = {0.1} / {size}; + Transfinite Curve {{5, 6}} = {0.02} / {size}; + Transfinite Curve {{-2, 4, 7}}={0.02}/{size} Using Bump 1/{bl_ratio}; + Transfinite Surface {{1, 2}} Right; + Mesh.MeshSizeExtendFromBoundary = 0; + Mesh.MeshSizeFromPoints = 0; + Mesh.MeshSizeFromCurvature = 0; + Mesh.Algorithm = 5; + Mesh.OptimizeNetgen = 1; + Mesh.Smoothing = 0; + """ + my_string = my_string + transfinite_string + else: + mesh_tail = \ + f""" + // Create distance field from surfaces, fluid sides + Field[1] = Distance; + Field[1].SurfacesList = {{1, 20, 28, 29}}; + Field[1].Sampling = 1000; + + //Create threshold field that varrries element size near boundaries + Field[2] = Threshold; + Field[2].InField = 1; + Field[2].SizeMin = {size} / {bl_ratio}; + Field[2].SizeMax = {size}; + Field[2].DistMin = 0.0002; + Field[2].DistMax = 0.010; + Field[2].StopAtDistMax = 1; + + // background mesh size + Field[3] = Box; + Field[3].XMin = 0.; + Field[3].XMax = 2.0; + Field[3].YMin = -1.0; + Field[3].YMax = 1.0; + Field[3].ZMin = -1.0; + Field[3].ZMax = 1.0; + Field[3].VIn = {size}; + + // Create distance field from surfaces, fluid/wall interface + Field[4] = Distance; + Field[4].SurfacesList = {{24}}; + Field[4].Sampling = 1000; + + //Create threshold field that varrries element size near boundaries + Field[5] = Threshold; + Field[5].InField = 4; + Field[5].SizeMin = {size} / {interface_ratio}; + Field[5].SizeMax = {size}; + Field[5].DistMin = 0.0002; + Field[5].DistMax = 0.010; + Field[5].StopAtDistMax = 1; + + // take the minimum of all defined meshing fields + Field[100] = Min; + Field[100].FieldsList = {{2, 3, 5}}; + Background Field = 100; + + Mesh.MeshSizeExtendFromBoundary = 0; + Mesh.MeshSizeFromPoints = 0; + Mesh.MeshSizeFromCurvature = 0; + + Mesh.Algorithm = 5; + Mesh.Algorithm3D = 10; + Mesh.OptimizeNetgen = 1; + Mesh.Smoothing = 100; + """ + my_string = my_string + mesh_tail + + #print(my_string) + from functools import partial + generate_mesh = partial(generate_gmsh, ScriptSource(my_string, "geo"), + force_ambient_dim=dim, dimensions=dim, target_unit="M", + return_tag_to_elements_map=True) + + return generate_mesh + + +def get_mesh_data(dim, mesh_angle, mesh_size, bl_ratio, interface_ratio, + transfinite, use_gmsh=True, periodic=False): + + # from meshmode.mesh.io import read_gmsh + # mesh, tag_to_elements = read_gmsh( + # mesh_filename, force_ambient_dim=dim, + # return_tag_to_elements_map=True) + if dim == 2: + mesh, tag_to_elements = get_flash_mesh_2d(angle=mesh_angle, + size=mesh_size, + bl_ratio=bl_ratio, + interface_ratio=interface_ratio, + transfinite=transfinite)() + else: + mesh, tag_to_elements = get_flash_mesh_3d(angle=mesh_angle, + size=mesh_size, + bl_ratio=bl_ratio, + interface_ratio=interface_ratio, + transfinite=transfinite)() + + volume_to_tags = { + "fluid": ["fluid"], + "wall": ["wall"]} + + # apply periodicity + if periodic: + from meshmode.mesh.processing import ( + glue_mesh_boundaries, BoundaryPairMapping) + + from meshmode import AffineMap + bdry_pair_mappings_and_tols = [] + offset = [0., 0.02] + bdry_pair_mappings_and_tols.append(( + BoundaryPairMapping( + "fluid_wall_bottom", + "fluid_wall_top", + AffineMap(offset=offset)), + 1e-12)) + + bdry_pair_mappings_and_tols.append(( + BoundaryPairMapping( + "solid_wall_bottom", + "solid_wall_top", + AffineMap(offset=offset)), + 1e-12)) + + mesh = glue_mesh_boundaries(mesh, bdry_pair_mappings_and_tols) + + return mesh, tag_to_elements, volume_to_tags + + +class InitSponge: + r"""Solution initializer for flow in the ACT-II facility + + This initializer creates a physics-consistent flow solution + given the top and bottom geometry profiles and an EOS using isentropic + flow relations. + + The flow is initialized from the inlet stagnations pressure, P0, and + stagnation temperature T0. + + geometry locations are linearly interpolated between given data points + + .. automethod:: __init__ + .. automethod:: __call__ + """ + def __init__(self, *, x0, thickness, amplitude): + r"""Initialize the sponge parameters. + + Parameters + ---------- + x0: float + sponge starting x location + thickness: float + sponge extent + amplitude: float + sponge strength modifier + """ + + self._x0 = x0 + self._thickness = thickness + self._amplitude = amplitude + + def __call__(self, x_vec, *, time=0.0): + """Create the sponge intensity at locations *x_vec*. + + Parameters + ---------- + x_vec: numpy.ndarray + Coordinates at which solution is desired + time: float + Time at which solution is desired. The strength is (optionally) + dependent on time + """ + xpos = x_vec[0] + actx = xpos.array_context + zeros = 0*xpos + x0 = zeros + self._x0 + + return self._amplitude * actx.np.where( + actx.np.greater(xpos, x0), + (zeros + ((xpos - self._x0)/self._thickness) + * ((xpos - self._x0)/self._thickness)), + zeros + 0.0 + ) + + +def getIsentropicPressure(mach, P0, gamma): + pressure = (1. + (gamma - 1.)*0.5*mach**2) + pressure = P0*pressure**(-gamma / (gamma - 1.)) + return pressure + + +def getIsentropicTemperature(mach, T0, gamma): + temperature = (1. + (gamma - 1.)*0.5*mach**2) + temperature = T0/temperature + return temperature + + +def smooth_step(actx, x, epsilon=1e-12): + # return actx.np.tanh(x) + # return actx.np.where( + # actx.np.greater(x, 0), + # actx.np.tanh(x)**3, + # 0*x) + return ( + actx.np.greater(x, 0) * actx.np.less(x, 1) * (1 - actx.np.cos(np.pi*x))/2 + + actx.np.greater(x, 1)) + + +class InitShock: + r"""Solution initializer for flow in the ACT-II facility + + This initializer creates a physics-consistent flow solution + given the top and bottom geometry profiles and an EOS using isentropic + flow relations. + + The flow is initialized from the inlet stagnations pressure, P0, and + stagnation temperature T0. + + geometry locations are linearly interpolated between given data points + + .. automethod:: __init__ + .. automethod:: __call__ + """ + + def __init__( + self, *, dim=2, nspecies=0, + P0, T0, temp_wall, temp_sigma, vel_sigma, gamma_guess, + mass_frac=None + ): + r"""Initialize mixture parameters. + + Parameters + ---------- + dim: int + specifies the number of dimensions for the solution + P0: float + stagnation pressure + T0: float + stagnation temperature + gamma_guess: float + guesstimate for gamma + temp_wall: float + wall temperature + temp_sigma: float + near-wall temperature relaxation parameter + vel_sigma: float + near-wall velocity relaxation parameter + """ + + if mass_frac is None: + if nspecies > 0: + mass_frac = np.zeros(shape=(nspecies,)) + + self._dim = dim + self._nspecies = nspecies + self._P0 = P0 + self._T0 = T0 + self._temp_wall = temp_wall + self._temp_sigma = temp_sigma + self._vel_sigma = vel_sigma + self._gamma_guess = gamma_guess + self._mass_frac = mass_frac + + def __call__(self, dcoll, x_vec, eos, *, time=0.0): + """Create the solution state at locations *x_vec*. + + Parameters + ---------- + x_vec: numpy.ndarray + Coordinates at which solution is desired + eos: + Mixture-compatible equation-of-state object must provide + these functions: + `eos.get_density` + `eos.get_internal_energy` + time: float + Time at which solution is desired. The location is (optionally) + dependent on time + """ + if x_vec.shape != (self._dim,): + raise ValueError(f"Position vector has unexpected dimensionality," + f" expected {self._dim}.") + + xpos = x_vec[0] + ypos = x_vec[1] + if self._dim == 3: + zpos = x_vec[2] + ytop = 0*x_vec[0] + actx = xpos.array_context + zeros = 0*xpos + ones = zeros + 1.0 + + mach = zeros + ytop = zeros + ybottom = zeros + # theta = zeros + gamma = self._gamma_guess + + pressure = getIsentropicPressure( + mach=mach, + P0=self._P0, + gamma=gamma + ) + temperature = getIsentropicTemperature( + mach=mach, + T0=self._T0, + gamma=gamma + ) + + # save the unsmoothed temerature, so we can use it with the injector init + # unsmoothed_temperature = temperature + + # modify the temperature in the near wall region to match the + # isothermal boundaries + sigma = self._temp_sigma + # wall_temperature = self._temp_wall + smoothing_top = smooth_step(actx, -sigma*(ypos-ytop)) + smoothing_bottom = smooth_step( + actx, sigma*actx.np.abs(ypos-ybottom)) + smoothing_fore = ones + smoothing_aft = ones + z0 = 0. + z1 = 0.035 + if self._dim == 3: + smoothing_fore = smooth_step(actx, sigma*(zpos-z0)) + smoothing_aft = smooth_step(actx, -sigma*(zpos-z1)) + + # smooth_temperature = (wall_temperature + + # (temperature - wall_temperature)*smoothing_top*smoothing_bottom * + # smoothing_fore*smoothing_aft) + + y = ones*self._mass_frac + mass = eos.get_density(pressure=pressure, temperature=temperature, + species_mass_fractions=y) + energy = mass*eos.get_internal_energy(temperature=temperature, + species_mass_fractions=y) + + velocity = np.zeros(self._dim, dtype=object) + mom = mass*velocity + + from mirgecom.fluid import make_conserved + cv = make_conserved(dim=self._dim, mass=mass, momentum=mom, energy=energy, + species_mass=mass*y) + velocity[0] = mach*eos.sound_speed(cv, temperature) + + # modify the velocity in the near-wall region to have a smooth profile + # this approximates the BL velocity profile + sigma = self._vel_sigma + smoothing_top = smooth_step(actx, -sigma*(ypos-ytop)) + smoothing_bottom = smooth_step(actx, sigma*(actx.np.abs(ypos-ybottom))) + smoothing_fore = ones + smoothing_aft = ones + if self._dim == 3: + smoothing_fore = smooth_step(actx, sigma*(zpos-z0)) + smoothing_aft = smooth_step(actx, -sigma*(zpos-z1)) + velocity[0] = (velocity[0]*smoothing_top*smoothing_bottom + * smoothing_fore*smoothing_aft) + + mom = mass*velocity + energy = (energy + np.dot(mom, mom)/(2.0*mass)) + return make_conserved( + dim=self._dim, + mass=mass, + momentum=mom, + energy=energy, + species_mass=mass*y + ) + + +def get_boundaries(dcoll, actx, dd_vol_fluid, dd_vol_wall, noslip, adiabatic, + periodic, temp_wall, gas_model, quadrature_tag, + target_fluid_state): + + from mirgecom.boundary import ( + PrescribedFluidBoundary, + IsothermalWallBoundary, + # SymmetryBoundary, + AdiabaticSlipBoundary, + AdiabaticNoslipWallBoundary, + DummyBoundary + ) + from mirgecom.diffusion import ( + DirichletDiffusionBoundary + ) + from mirgecom.simutil import force_evaluation + from mirgecom.gas_model import project_fluid_state + + if noslip: + if adiabatic: + fluid_wall = AdiabaticNoslipWallBoundary() + else: + fluid_wall = IsothermalWallBoundary(temp_wall) + + else: + # new impl, following Mengaldo with modifications for slip vs no slip + # set the flux directly, instead using viscous numerical flux func + # fluid_wall = AdiabaticSlipWallBoundary2() + + # implementation from mirgecom + # should be same as AdiabaticSlipBoundary2 + fluid_wall = AdiabaticSlipBoundary() + + # new impl, following Mengaldo with modifications for slip vs no slip + # local version + # fluid_wall = AdiabaticSlipWallBoundary() + + # Tulio's symmetry boundary + # fluid_wall = SymmetryBoundary(dim=dim) + + wall_farfield = DirichletDiffusionBoundary(temp_wall) + + # use dummy boundaries to setup the smoothness state for the target + target_boundaries = { + dd_vol_fluid.trace("fluid_inflow").domain_tag: DummyBoundary(), + # dd_vol_fluid.trace("fluid_wall").domain_tag: IsothermalWallBoundary() + dd_vol_fluid.trace("fluid_wall").domain_tag: fluid_wall + } + + def get_target_state_on_boundary(btag): + return project_fluid_state( + dcoll, dd_vol_fluid, + dd_vol_fluid.trace(btag).with_discr_tag(quadrature_tag), + target_fluid_state, gas_model + ) + + flow_ref_state = \ + get_target_state_on_boundary("fluid_inflow") + + flow_ref_state = force_evaluation(actx, flow_ref_state) + + def _target_flow_state_func(**kwargs): + return flow_ref_state + + flow_boundary = PrescribedFluidBoundary( + boundary_state_func=_target_flow_state_func) + + if periodic: + fluid_boundaries = { + dd_vol_fluid.trace("fluid_inflow").domain_tag: flow_boundary, + } + + wall_boundaries = { + dd_vol_wall.trace("solid_wall_end").domain_tag: wall_farfield + } + else: + fluid_boundaries = { + dd_vol_fluid.trace("fluid_inflow").domain_tag: flow_boundary, + dd_vol_fluid.trace("fluid_wall").domain_tag: fluid_wall + } + + wall_boundaries = { + dd_vol_wall.trace("wall_farfield").domain_tag: wall_farfield + } + return fluid_boundaries, wall_boundaries, target_boundaries diff --git a/y3prediction/prediction.py b/y3prediction/prediction.py index 0a71e70..73cbba8 100644 --- a/y3prediction/prediction.py +++ b/y3prediction/prediction.py @@ -25,6 +25,7 @@ """ import logging import sys +import gc import numpy as np import pyopencl as cl import numpy.linalg as la # noqa @@ -71,13 +72,6 @@ from mirgecom.fluid import make_conserved from mirgecom.limiter import bound_preserving_limiter from mirgecom.steppers import advance_state -from mirgecom.boundary import ( - PrescribedFluidBoundary, - IsothermalWallBoundary, - AdiabaticSlipBoundary, - AdiabaticNoslipWallBoundary, - DummyBoundary -) from mirgecom.diffusion import ( diffusion_operator, DirichletDiffusionBoundary, @@ -87,6 +81,7 @@ from mirgecom.eos import IdealSingleGas, PyrometheusMixture from mirgecom.transport import (SimpleTransport, PowerLawTransport, + MixtureAveragedTransport, ArtificialViscosityTransportDiv, ArtificialViscosityTransportDiv2) from mirgecom.gas_model import GasModel, make_fluid_state @@ -271,13 +266,14 @@ def main(ctx_factory=cl.create_some_context, from mirgecom.io import read_and_distribute_yaml_data input_data = read_and_distribute_yaml_data(comm, user_input_file) + # Let the user specify a name to customize initialization + init_name = configurate("init_name", input_data, "ACTII") + # i/o frequencies nviz = configurate("nviz", input_data, 500) nrestart = configurate("nrestart", input_data, 5000) nhealth = configurate("nhealth", input_data, 1) nstatus = configurate("nstatus", input_data, 1) - - # garbage collection frequency ngarbage = configurate("ngarbage", input_data, 10) # verbosity for what gets written to viz dumps, increase for more stuff @@ -315,10 +311,10 @@ def main(ctx_factory=cl.create_some_context, alpha_sc = configurate("alpha_sc", input_data, 0.3) kappa_sc = configurate("kappa_sc", input_data, 0.5) s0_sc = configurate("s0_sc", input_data, -5.0) - av2_mu0 = configurate("av_mu0", input_data, 0.1) + av2_mu0 = configurate("av2_mu0", input_data, 0.1) av2_beta0 = configurate("av2_beta0", input_data, 6.0) av2_kappa0 = configurate("av2_kappa0", input_data, 1.0) - av2_prandtl0 = configurate("av2_prandtl0", input_data, 0.9) + av2_prandtl0 = configurate("av_prandtl0", input_data, 0.9) av2_mu_s0 = configurate("av2_mu_s0", input_data, 0.) av2_kappa_s0 = configurate("av2_kappa_s0", input_data, 0.) av2_beta_s0 = configurate("av2_beta_s0", input_data, 0.01) @@ -339,6 +335,67 @@ def main(ctx_factory=cl.create_some_context, adiabatic = configurate("adiabatic", input_data, False) use_1d_part = configurate("use_1d_part", input_data, True) + # - Flash1D-specific configuration parameters + # -- mesh generation + mesh_angle = configurate("mesh_angle", input_data, 0.) + mesh_size = configurate("mesh_size", input_data, 1e-3) + bl_ratio = configurate("bl_ratio", input_data, 3.) + interface_ratio = configurate("interface_ratio", input_data, 2.) + transfinite = configurate("transfinite", input_data, False) + use_gmsh = configurate("use_gmsh", input_data, True) + periodic = configurate("periodic", input_data, False) + mesh_theta = mesh_angle/180.*np.pi/2. + mesh_normal = np.zeros(shape=(dim,)) + mesh_normal[0] = np.cos(mesh_theta) + mesh_normal[1] = np.sin(mesh_theta) + mesh_normal = mesh_normal / np.sqrt(np.dot(mesh_normal, mesh_normal)) + mesh_rotation = np.matrix([[mesh_normal[0], + mesh_normal[1]], + [-mesh_normal[1], + mesh_normal[0]]]) + # -- initialization parameters + # --- common + flash_sigma = configurate("sigma", input_data, .001) + flash_mach = configurate("mach", input_data, 2.0) + flash_press = configurate("pressure_background", input_data, 101325.) + flash_temp = configurate("temperature_background", input_data, 300.) + init_flame = configurate("init_flame", input_data, False) + init_shock = configurate("init_shock", input_data, False) + + # --- shock init + shock_angle = configurate("shock_angle", input_data, 0.) + shock_normal = np.zeros(shape=(dim,)) + shock_theta = shock_angle/180.*np.pi/2. + shock_normal[0] = np.cos(shock_theta) + shock_normal[1] = np.sin(shock_theta) + shock_normal = shock_normal / np.sqrt(np.dot(shock_normal, shock_normal)) + shock_location = np.zeros(shape=(dim,)) + shock_location[0] = configurate("shock_location_x", input_data, .05) + shock_location[1] = configurate("shock_location_y", input_data, 0.) + if not init_shock: + shock_location = None + + # --- flame init + flame_angle = configurate("flame_angle", input_data, 0.) + flame_normal = np.zeros(shape=(dim,)) + flame_theta = flame_angle/180.*np.pi/2. + flame_normal[0] = np.cos(flame_theta) + flame_normal[1] = np.sin(flame_theta) + flame_location = np.zeros(shape=(dim,)) + flame_location[0] = configurate("flame_location_x", input_data, .05) + flame_location[1] = configurate("flame_location_y", input_data, 0.) + temperature_burned = configurate("temperature_burned", input_data, 1500.) + if not init_flame: + flame_location = None + + if np.abs(mesh_angle) > 0.: + if init_shock: + shock_location = mesh_rotation.dot(shock_location) + shock_normal = mesh_rotation.dot(shock_normal) + if init_flame: + flame_location = mesh_rotation.dot(flame_location) + flame_normal = mesh_rotation.dot(flame_normal) + # material properties and models options gas_mat_prop = configurate("gas_mat_prop", input_data, 0) spec_diff = configurate("spec_diff", input_data, 1.e-4) @@ -376,9 +433,7 @@ def main(ctx_factory=cl.create_some_context, use_av = configurate("use_av", input_data, 0) # species limiter - # 0 - none - # 1 - limit in on call to make_fluid_state - use_species_limiter = configurate("use_species_limiter", input_data, 0) + use_species_limiter = configurate("use_species_limiter", input_data, False) # Filtering is implemented according to HW Sec. 5.3 # The modal response function is e^-(alpha * eta ^ 2s), where @@ -423,21 +478,21 @@ def main(ctx_factory=cl.create_some_context, total_temp_inflow = configurate("total_temp_inflow", input_data, 2076.43) # injection flow properties - total_pres_inj = configurate("total_pres_inj", input_data, 50400) - total_temp_inj = configurate("total_temp_inj", input_data, 300) + total_pres_inj = configurate("total_pres_inj", input_data, 50400.) + total_temp_inj = configurate("total_temp_inj", input_data, 300.) mach_inj = configurate("mach_inj", input_data, 1.0) # parameters to adjust the shape of the initialization - vel_sigma = configurate("vel_sigma", input_data, 1000) - temp_sigma = configurate("temp_sigma", input_data, 1250) + vel_sigma = configurate("vel_sigma", input_data, 1000.) + temp_sigma = configurate("temp_sigma", input_data, 1250.) # adjusted to match the mass flow rate - vel_sigma_inj = configurate("vel_sigma_inj", input_data, 5000) - temp_sigma_inj = configurate("temp_sigma_inj", input_data, 5000) + vel_sigma_inj = configurate("vel_sigma_inj", input_data, 5000.) + temp_sigma_inj = configurate("temp_sigma_inj", input_data, 5000.) temp_wall = 300 # wall stuff - wall_penalty_amount = configurate("wall_penalty_amount", input_data, 0) - wall_time_scale = configurate("wall_time_scale", input_data, 1) + wall_penalty_amount = configurate("wall_penalty_amount", input_data, 0.) + wall_time_scale = configurate("wall_time_scale", input_data, 1.) wall_material = configurate("wall_material", input_data, 0) # use fluid average diffusivity by default @@ -445,15 +500,15 @@ def main(ctx_factory=cl.create_some_context, # Averaging from https://www.azom.com/article.aspx?ArticleID=1630 # for graphite - wall_insert_rho = configurate("wall_insert_rho", input_data, 1625) - wall_insert_cp = configurate("wall_insert_cp", input_data, 770) + wall_insert_rho = configurate("wall_insert_rho", input_data, 1625.) + wall_insert_cp = configurate("wall_insert_cp", input_data, 770.) wall_insert_kappa = configurate("wall_insert_kappa", input_data, 247.5) # Averaging from http://www.matweb.com/search/datasheet.aspx?bassnum=MS0001 # for steel wall_surround_rho = configurate("wall_surround_rho", input_data, 7.9e3) - wall_surround_cp = configurate("wall_surround_cp", input_data, 470) - wall_surround_kappa = configurate("wall_surround_kappa", input_data, 48) + wall_surround_cp = configurate("wall_surround_cp", input_data, 470.) + wall_surround_kappa = configurate("wall_surround_kappa", input_data, 48.) # initialize the ignition spark spark_init_loc_z = 0.035/2. @@ -483,7 +538,6 @@ def main(ctx_factory=cl.create_some_context, print(f"Shock capturing parameters: alpha {alpha_sc}, " f"s0 {s0_sc}, kappa {kappa_sc}") - # use_av=1 specific parameters # flow stagnation temperature static_temp = 2076.43 # steepness of the smoothed function @@ -514,12 +568,12 @@ def main(ctx_factory=cl.create_some_context, print("Artificial viscosity using modified transport properties") print("\t mu, beta, kappa") # MJA update this - print(f"Shock capturing parameters:" - f"\tav_mu {av2_mu0}" - f"\tav_beta {av2_beta0}" - f"\tav_kappa {av2_kappa0}" - f"\tav_prantdl {av2_prandtl0}" - f"stagnation temperature {static_temp}") + print(f"Shock capturing parameters:\n" + f"\tav_mu {av2_mu0}\n" + f"\tav_beta {av2_beta0}\n" + f"\tav_kappa {av2_kappa0}\n" + f"\tav_prantdl {av2_prandtl0}\n" + f"\tstagnation temperature {static_temp}\n") else: error_message = "Unknown artifical viscosity model {}".format(use_av) raise RuntimeError(error_message) @@ -529,6 +583,11 @@ def main(ctx_factory=cl.create_some_context, print(f"\tnrestart = {nrestart}") print(f"\tnhealth = {nhealth}") print(f"\tnstatus = {nstatus}") + if ngarbage >= 0: + print(f"\tSyncd garbage collection every {ngarbage} steps.") + # gc.disable() + else: + print("\tUsing Python automatic garbage collection.") if constant_cfl == 1: print(f"\tConstant cfl mode, current_cfl = {current_cfl}") else: @@ -538,9 +597,9 @@ def main(ctx_factory=cl.create_some_context, print(f"\tdimension = {dim}") print(f"\tTime integration {integrator}") if noslip: - print("Fluid wall boundary conditions are noslip for veloctiy") + print("Fluid wall boundary conditions are noslip for velocity") else: - print("Fluid wall boundary conditions are slip for veloctiy") + print("Fluid wall boundary conditions are slip for velocity") if adiabatic: print("Fluid wall boundary conditions are adiabatic for temperature") else: @@ -575,34 +634,36 @@ def main(ctx_factory=cl.create_some_context, if adiabatic: temp_sigma = 0. """ - - if rank == 0: - print("\n#### Simluation setup data: ####") - print(f"\ttotal_pres_injection = {total_pres_inj}") - print(f"\ttotal_temp_injection = {total_temp_inj}") - print(f"\tvel_sigma = {vel_sigma}") - print(f"\ttemp_sigma = {temp_sigma}") - print(f"\tvel_sigma_injection = {vel_sigma_inj}") - print(f"\ttemp_sigma_injection = {temp_sigma_inj}") - print("#### Simluation setup data: ####") + if init_name == "ACTII": + if rank == 0: + print("\n#### Simluation setup data: ####") + print(f"\ttotal_pres_injection = {total_pres_inj}") + print(f"\ttotal_temp_injection = {total_temp_inj}") + print(f"\tvel_sigma = {vel_sigma}") + print(f"\ttemp_sigma = {temp_sigma}") + print(f"\tvel_sigma_injection = {vel_sigma_inj}") + print(f"\ttemp_sigma_injection = {temp_sigma_inj}") + print("#### Simluation setup data: ####") spark_center = np.zeros(shape=(dim,)) spark_center[0] = spark_init_loc_x spark_center[1] = spark_init_loc_y if dim == 3: spark_center[2] = spark_init_loc_z - if rank == 0 and use_ignition > 0: - print("\n#### Ignition control parameters ####") - print(f"spark center ({spark_center[0]},{spark_center[1]})") - print(f"spark FWHM {spark_diameter}") - print(f"spark strength {spark_strength}") - print(f"ignition time {spark_init_time}") - print(f"ignition duration {spark_duration}") - if use_ignition == 1: - print("spark ignition") - elif use_ignition == 2: - print("heat source ignition") - print("#### Ignition control parameters ####\n") + + if init_name == "ACTII": + if rank == 0 and use_ignition > 0: + print("\n#### Ignition control parameters ####") + print(f"spark center ({spark_center[0]},{spark_center[1]})") + print(f"spark FWHM {spark_diameter}") + print(f"spark strength {spark_strength}") + print(f"ignition time {spark_init_time}") + print(f"ignition duration {spark_duration}") + if use_ignition == 1: + print("spark ignition") + elif use_ignition == 2: + print("heat source ignition") + print("#### Ignition control parameters ####\n") def _compiled_stepper_wrapper(state, t, dt, rhs): return compiled_lsrk45_step(actx, state, t, dt, rhs) @@ -691,7 +752,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # don't allow limiting on flows without species if nspecies == 0: - use_species_limiter = 0 + use_species_limiter = False use_injection = False # Turn off combustion unless EOS supports it @@ -728,7 +789,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): elif eos_type == 1: print("\tPyrometheus EOS") - if use_species_limiter == 1: + if use_species_limiter: print("\nSpecies mass fractions limited to [0:1]") transport_alpha = 0.6 @@ -776,9 +837,9 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): print("\tWall oxidizer transport disabled") if use_wall_mass: - print("\t Wall mass loss enabled") + print("\tWall mass loss enabled") else: - print("\t Wall mass loss disabled") + print("\tWall mass loss disabled") print(f"\tWall density = {wall_insert_rho}") print(f"\tWall cp = {wall_insert_cp}") @@ -790,6 +851,13 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): print(f"\tWall penalty = {wall_penalty_amount}") print("#### Simluation material properties: ####") + if eos_type == 1 or transport_type == 2: + from mirgecom.thermochemistry import get_pyrometheus_wrapper_class + chem_source_tol = 1.e-10 + pyro_mech = get_pyrometheus_wrapper_class( + pyro_class=Thermochemistry, temperature_niter=pyro_temp_iter, + zero_level=chem_source_tol)(actx.np) + spec_diffusivity = spec_diff * np.ones(nspecies) if transport_type == 0: physical_transport_model = SimpleTransport( @@ -800,6 +868,8 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): alpha=transport_alpha, beta=transport_beta, sigma=transport_sigma, n=transport_n, species_diffusivity=spec_diffusivity) + if transport_type == 2: + physical_transport_model = MixtureAveragedTransport(pyro_mech) transport_model = physical_transport_model if use_av == 1: @@ -829,22 +899,19 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): inlet_area_ratio = inlet_height/throat_height outlet_area_ratio = outlet_height/throat_height - chem_source_tol = 1.e-10 # make the eos if eos_type == 0: eos = IdealSingleGas(gamma=gamma, gas_const=r) species_names = ["air", "fuel", "inert"] else: - from mirgecom.thermochemistry import get_pyrometheus_wrapper_class - pyro_mech = get_pyrometheus_wrapper_class( - pyro_class=Thermochemistry, temperature_niter=pyro_temp_iter, - zero_level=chem_source_tol)(actx.np) eos = PyrometheusMixture(pyro_mech, temperature_guess=init_temperature) species_names = pyro_mech.species_names gas_model = GasModel(eos=eos, transport=transport_model) - # initialize eos and species mass fractions + # initialize species mass fractions + inj_ymin = -0.0243245 + inj_ymax = -0.0227345 y = np.zeros(nspecies) y_fuel = np.zeros(nspecies) if nspecies == 3: @@ -869,6 +936,23 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): y_fuel[i_c2h4] = mf_c2h4 y_fuel[i_h2] = mf_h2 + # + # stagnation tempertuare 2076.43 K + # stagnation pressure 2.745e5 Pa + # + # isentropic expansion based on the area ratios between the inlet (r=54e-3m) and + # the throat (r=3.167e-3) + # + vel_inflow = np.zeros(shape=(dim,)) + vel_outflow = np.zeros(shape=(dim,)) + vel_injection = np.zeros(shape=(dim,)) + + throat_height = 3.61909e-3 + inlet_height = 54.129e-3 + outlet_height = 28.54986e-3 + inlet_area_ratio = inlet_height/throat_height + outlet_area_ratio = outlet_height/throat_height + inlet_mach = getMachFromAreaRatio(area_ratio=inlet_area_ratio, gamma=gamma, mach_guess=0.01) @@ -919,15 +1003,16 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): vel_inflow[0] = inlet_mach*sos - if rank == 0: - print("#### Simluation initialization data: ####") - print(f"\tinlet Mach number {inlet_mach}") - print(f"\tinlet gamma {inlet_gamma}") - print(f"\tinlet temperature {temp_inflow}") - print(f"\tinlet pressure {pres_inflow}") - print(f"\tinlet rho {rho_inflow}") - print(f"\tinlet velocity {vel_inflow[0]}") - #print(f"final inlet pressure {pres_inflow_final}") + if init_name == "ACTII": + if rank == 0: + print("#### Simluation initialization data: ####") + print(f"\tinlet Mach number {inlet_mach}") + print(f"\tinlet gamma {inlet_gamma}") + print(f"\tinlet temperature {temp_inflow}") + print(f"\tinlet pressure {pres_inflow}") + print(f"\tinlet rho {rho_inflow}") + print(f"\tinlet velocity {vel_inflow[0]}") + #print(f"final inlet pressure {pres_inflow_final}") outlet_mach = getMachFromAreaRatio(area_ratio=outlet_area_ratio, gamma=gamma, @@ -977,14 +1062,15 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): vel_outflow[0] = outlet_mach*math.sqrt(gamma*pres_outflow/rho_outflow) - if rank == 0: - print("\t********") - print(f"\toutlet Mach number {outlet_mach}") - print(f"\toutlet gamma {outlet_gamma}") - print(f"\toutlet temperature {temp_outflow}") - print(f"\toutlet pressure {pres_outflow}") - print(f"\toutlet rho {rho_outflow}") - print(f"\toutlet velocity {vel_outflow[0]}") + if init_name == "ACTII": + if rank == 0: + print("\t********") + print(f"\toutlet Mach number {outlet_mach}") + print(f"\toutlet gamma {outlet_gamma}") + print(f"\toutlet temperature {temp_outflow}") + print(f"\toutlet pressure {pres_outflow}") + print(f"\toutlet rho {rho_outflow}") + print(f"\toutlet velocity {vel_outflow[0]}") gamma_injection = gamma if nspecies > 0: @@ -1042,49 +1128,136 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): sos = math.sqrt(gamma_injection*pres_injection/rho_injection) vel_injection[0] = -mach_inj*sos - - if rank == 0: - print("\t********") - print(f"\tinjector Mach number {mach_inj}") - print(f"\tinjector gamma {gamma_injection}") - print(f"\tinjector temperature {temp_injection}") - print(f"\tinjector pressure {pres_injection}") - print(f"\tinjector rho {rho_injection}") - print(f"\tinjector velocity {vel_injection[0]}") - print("#### Simluation initialization data: ####\n") + if init_name == "ACTII": + if rank == 0: + print("\t********") + print(f"\tinjector Mach number {mach_inj}") + print(f"\tinjector gamma {gamma_injection}") + print(f"\tinjector temperature {temp_injection}") + print(f"\tinjector pressure {pres_injection}") + print(f"\tinjector rho {rho_injection}") + print(f"\tinjector velocity {vel_injection[0]}") + print("#### Simluation initialization data: ####\n") else: - if rank == 0: - print("\t********") - print("\tnspecies=0, injection disabled") + if init_name == "ACTII": + if rank == 0: + print("\t********") + print("\tnspecies=0, injection disabled") # read geometry files geometry_bottom = None geometry_top = None - if rank == 0: - from numpy import loadtxt - geometry_bottom = loadtxt("data/nozzleBottom.dat", - comments="#", unpack=False) - geometry_top = loadtxt("data/nozzleTop.dat", - comments="#", unpack=False) - geometry_bottom = comm.bcast(geometry_bottom, root=0) - geometry_top = comm.bcast(geometry_top, root=0) - inj_ymin = -0.0243245 - inj_ymax = -0.0227345 - bulk_init = InitACTII(dim=dim, - geom_top=geometry_top, geom_bottom=geometry_bottom, - P0=total_pres_inflow, T0=total_temp_inflow, - temp_wall=temp_wall, temp_sigma=temp_sigma, - vel_sigma=vel_sigma, nspecies=nspecies, - mass_frac=y, gamma_guess=inlet_gamma, - inj_gamma_guess=gamma_injection, - inj_pres=total_pres_inj, - inj_temp=total_temp_inj, - inj_vel=vel_injection, inj_mass_frac=y_fuel, - inj_temp_sigma=temp_sigma_inj, - inj_vel_sigma=vel_sigma_inj, - inj_ytop=inj_ymax, inj_ybottom=inj_ymin, - inj_mach=mach_inj, injection=use_injection) + if init_name == "ACTII": + if rank == 0: + print("\tInitializing ACTII domain.\n") + + from numpy import loadtxt + geometry_bottom = loadtxt("data/nozzleBottom.dat", + comments="#", unpack=False) + geometry_top = loadtxt("data/nozzleTop.dat", + comments="#", unpack=False) + + geometry_bottom = comm.bcast(geometry_bottom, root=0) + geometry_top = comm.bcast(geometry_top, root=0) + fluid_init = InitACTII(dim=dim, + geom_top=geometry_top, geom_bottom=geometry_bottom, + P0=total_pres_inflow, T0=total_temp_inflow, + temp_wall=temp_wall, temp_sigma=temp_sigma, + vel_sigma=vel_sigma, nspecies=nspecies, + mass_frac=y, gamma_guess=inlet_gamma, + inj_gamma_guess=gamma_injection, + inj_pres=total_pres_inj, + inj_temp=total_temp_inj, + inj_vel=vel_injection, inj_mass_frac=y_fuel, + inj_temp_sigma=temp_sigma_inj, + inj_vel_sigma=vel_sigma_inj, + inj_ytop=inj_ymax, inj_ybottom=inj_ymin, + inj_mach=mach_inj, injection=use_injection) + + elif init_name == "Flash1D": + print("\tInitializing flame and shock (Flash1D) domain.\n") + from y3prediction.flash_utils import Flash1D + rho_bkrnd = flash_press/r/flash_temp + c_bkrnd = math.sqrt(gamma*flash_press/rho_bkrnd) + + pressure_ratio = (2.*gamma*flash_mach**2-(gamma-1.))/(gamma+1.) + density_ratio = (gamma+1.)*flash_mach**2/((gamma-1.)*flash_mach**2+2.) + + rho1 = rho_bkrnd + pressure1 = flash_press + temperature1 = pressure1/rho1/r + rho2 = rho1*density_ratio + pressure2 = pressure1*pressure_ratio + temperature2 = pressure2/rho2/r + velocity2 = -flash_mach*c_bkrnd*(1/density_ratio-1) + temp_wall = temperature1 + + vel_left = velocity2*shock_normal + vel_right = np.zeros(shape=(dim,)) + vel_cross = np.zeros(shape=(dim,)) + vel_cross[1] = 0 + + pressure_shocked = pressure2 + pressure_unshocked = pressure1 + temperature_shocked = temperature2 + temperature_unshocked = temperature1 + vel_shocked = vel_left + vel_unshocked = vel_right + + if init_flame and nspecies == 7: + from y3prediction.flash_utils import setup_flame + flame_states = setup_flame(pressure_unburned=pressure_unshocked, + temperature_unburned=temperature_unshocked) + pressure_unburned, temperature_unburned, y_unburned, rho_unburned = \ + flame_states["unburned"] + pressure_burned, temperature_burned_c, y_burned, rho_burned = \ + flame_states["burned"] + else: + y_unburned = y_fuel + y_burned = y + pressure_burned = pressure_unshocked + pressure_unburned = pressure_unshocked + temperature_burned = temperature_unshocked + temperature_unburned = temperature_unshocked + + if rank == 0: + print("#### Flash1D initialization data: ####") + print(f"\tShock Mach number {flash_mach}") + print(f"\tgamma {gamma}") + print(f"\tambient temperature {temperature_unshocked}") + print(f"\tambient pressure {pressure_unshocked}") + print(f"\tambient rho {rho1}") + print(f"\tambient velocity {vel_unshocked}") + print(f"\tpost-shock temperature {temperature_shocked}") + print(f"\tpost-shock pressure {pressure_shocked}") + print(f"\tpost-shock rho {rho2}") + print(f"\tpost-shock velocity {vel_shocked}") + + fluid_init = Flash1D(dim=dim, + nspecies=nspecies, + shock_location=shock_location, + flame_location=flame_location, + shock_normal_dir=shock_normal, + flame_normal_dir=flame_normal, + sigma=flash_sigma, + pressure_shocked=pressure_shocked, + pressure_unshocked=pressure_unshocked, + temperature_shocked=temperature_shocked, + temperature_unshocked=temperature_unshocked, + velocity_shocked=vel_shocked, + velocity_unshocked=vel_unshocked, + velocity_cross=vel_cross, + species_mass_fractions_shocked=y, + species_mass_fractions_unshocked=y, + species_mass_fractions_burned=y_burned, + species_mass_fractions_unburned=y_unburned, + temperature_burned=temperature_burned, + temperature_unburned=temperature_unburned, + pressure_burned=pressure_burned, + pressure_unburned=pressure_unburned, + temp_wall=temp_wall, + vel_sigma=vel_sigma, temp_sigma=temp_sigma) viz_path = "viz_data/" vizname = viz_path + casename @@ -1110,19 +1283,23 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): assert restart_data["nparts"] == nparts assert restart_data["nspecies"] == nspecies - else: # generate the grid from scratch + else: # initialize grid data if rank == 0: - print(f"Reading mesh from {mesh_filename}") - - def get_mesh_data(): - from meshmode.mesh.io import read_gmsh - mesh, tag_to_elements = read_gmsh( - mesh_filename, force_ambient_dim=dim, - return_tag_to_elements_map=True) - volume_to_tags = { - "fluid": ["fluid"], - "wall": ["wall_insert", "wall_surround"]} - return mesh, tag_to_elements, volume_to_tags + print(f"Initializing grid data for initialization case: {init_name}") + if init_name == "ACTII": + from y3prediction.actii_y3 import get_mesh_data as get_init_mesh + get_mesh_data = \ + partial(get_init_mesh, mesh_filename=mesh_filename, dim=dim) + elif init_name == "Flash1D": + from y3prediction.flash_utils import get_mesh_data as get_init_mesh + get_mesh_data = \ + partial(get_init_mesh, dim=dim, mesh_angle=mesh_angle, + mesh_size=mesh_size, bl_ratio=bl_ratio, + interface_ratio=interface_ratio, + transfinite=transfinite, periodic=periodic, + use_gmsh=use_gmsh) + else: + raise ValueError(f"Unrecognized initializtion case: {init_name}.") def my_partitioner(mesh, tag_to_elements, num_ranks): from mirgecom.simutil import geometric_mesh_partitioner @@ -1181,17 +1358,30 @@ def my_partitioner(mesh, tag_to_elements, num_ranks): wall_vol_discr = dcoll.discr_from_dd(dd_vol_wall) wall_tag_to_elements = volume_to_local_mesh_data["wall"][1] - wall_insert_mask = mask_from_elements( - wall_vol_discr, actx, wall_tag_to_elements["wall_insert"]) - wall_surround_mask = mask_from_elements( - wall_vol_discr, actx, wall_tag_to_elements["wall_surround"]) - - wall_bnd = dd_vol_fluid.trace("isothermal_wall") - inflow_bnd = dd_vol_fluid.trace("inflow") - outflow_bnd = dd_vol_fluid.trace("outflow") - inj_bnd = dd_vol_fluid.trace("injection") - flow_bnd = dd_vol_fluid.trace("flow") - wall_ffld_bnd = dd_vol_wall.trace("wall_farfield") + + if init_name == "ACTII": + wall_insert_mask = mask_from_elements( + wall_vol_discr, actx, wall_tag_to_elements["wall_insert"]) + wall_surround_mask = mask_from_elements( + wall_vol_discr, actx, wall_tag_to_elements["wall_surround"]) + elif init_name == "Flash1D": + wall_insert_mask = mask_from_elements( + wall_vol_discr, actx, wall_tag_to_elements["wall"]) + wall_surround_mask = None + + # initialize the wall heat capacity + wall_cp = wall_insert_cp * wall_insert_mask + if wall_surround_mask is not None: + wall_cp = wall_cp + wall_surround_cp * wall_surround_mask + + if init_name == "ACTII": + wall_bnd = dd_vol_fluid.trace("isothermal_wall") + flow_bnd = dd_vol_fluid.trace("flow") + wall_ffld_bnd = dd_vol_wall.trace("wall_farfield") + elif init_name == "Flash1D": + wall_bnd = dd_vol_fluid.trace("fluid_wall") + flow_bnd = dd_vol_fluid.trace("fluid_inflow") + wall_ffld_bnd = dd_vol_wall.trace("wall_farfield") from grudge.dt_utils import characteristic_lengthscales char_length_fluid = force_evaluation(actx, @@ -1315,9 +1505,9 @@ def limit_fluid_state(cv, pressure, temperature, dd=dd_vol_fluid): if rhs_filter_cutoff < 0: rhs_filter_cutoff = int(rhs_filter_frac * order) - if soln_filter_cutoff >= order: + if soln_nfilter > 0 and soln_filter_cutoff >= order: raise ValueError("Invalid setting for solution filter (cutoff >= order).") - if rhs_filter_cutoff >= order: + if use_rhs_filter and rhs_filter_cutoff >= order: raise ValueError("Invalid setting for RHS filter (cutoff >= order).") from mirgecom.filter import ( @@ -1357,16 +1547,14 @@ def filter_wall_rhs(rhs): logger.info(f" - filter cutoff = {rhs_filter_cutoff}") logger.info(f" - filter order = {rhs_filter_order}") - limiter_func = None - if use_species_limiter: - limiter_func = limit_fluid_state + limiter_func = limit_fluid_state if use_species_limiter else None ######################################## # Helper functions for building states # ######################################## - def _create_fluid_state(cv, temperature_seed, smoothness_mu=None, - smoothness_beta=None, smoothness_kappa=None): + def _create_fluid_state(cv, temperature_seed, smoothness_mu, + smoothness_beta, smoothness_kappa): return make_fluid_state(cv=cv, gas_model=gas_model, temperature_seed=temperature_seed, smoothness_mu=smoothness_mu, @@ -1442,7 +1630,7 @@ def compute_smoothness(cv, dv, grad_cv): compute_smoothness_compiled = actx.compile(compute_smoothness) # noqa def lmax(s): - b = 100 + b = 1000 return (s/np.pi*actx.np.arctan(b*s) + 0.5*s - 1/np.pi*actx.np.arctan(b) + 0.5) @@ -1575,7 +1763,6 @@ def grad_t_operator_fluid(fluid_state, time): ################################## # Set up flow initial conditions # ################################## - if restart_filename: if rank == 0: logger.info("Restarting soln.") @@ -1616,32 +1803,36 @@ def grad_t_operator_fluid(fluid_state, time): # Set the current state from time 0 if rank == 0: logger.info("Initializing soln.") - restart_cv = bulk_init( - dcoll=dcoll, x_vec=fluid_nodes, eos=eos, - time=0) + + if init_name == "ACTII": + restart_cv = fluid_init(dcoll=dcoll, x_vec=fluid_nodes, eos=eos, time=0) + else: + restart_cv = fluid_init(x_vec=fluid_nodes, eos=eos, time=0) restart_cv = force_evaluation(actx, restart_cv) temperature_seed = actx.zeros_like(restart_cv.mass) + init_temperature temperature_seed = force_evaluation(actx, temperature_seed) # create a fluid state so we can compute grad_t and grad_cv - restart_fluid_state = create_fluid_state(cv=restart_cv, - temperature_seed=temperature_seed) restart_av_smu = actx.zeros_like(restart_cv.mass) restart_av_sbeta = actx.zeros_like(restart_cv.mass) restart_av_skappa = actx.zeros_like(restart_cv.mass) + restart_fluid_state = create_fluid_state(cv=restart_cv, + temperature_seed=temperature_seed, + smoothness_mu=restart_av_smu, + smoothness_beta=restart_av_sbeta, + smoothness_kappa=restart_av_skappa) # Ideally we would compute the smoothness variables here, # but we need the boundary conditions (and hence the target state) first, # so we defer until after those are setup # initialize the wall - wall_mass = ( - wall_insert_rho * wall_insert_mask - + wall_surround_rho * wall_surround_mask) - wall_cp = ( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask) + wall_mass = wall_insert_rho * wall_insert_mask + + if wall_surround_mask is not None: + wall_mass = wall_mass + wall_surround_rho * wall_surround_mask + restart_wv = WallVars( mass=wall_mass, energy=wall_mass * wall_cp * temp_wall, @@ -1718,124 +1909,47 @@ def grad_t_operator_target(fluid_state, time): # use dummy boundaries to update the smoothness state for the target if use_av > 0: - if use_injection: - target_boundaries = { - flow_bnd.domain_tag: # pylint: disable=no-member - DummyBoundary(), - wall_bnd.domain_tag: # pylint: disable=no-member - IsothermalWallBoundary() - } - else: - target_boundaries = { - inflow_bnd.domain_tag: # pylint: disable=no-member - DummyBoundary(), - outflow_bnd.domain_tag: # pylint: disable=no-member - DummyBoundary(), - inj_bnd.domain_tag: # pylint: disable=no-member - IsothermalWallBoundary(), - wall_bnd.domain_tag: # pylint: disable=no-member - IsothermalWallBoundary() - } - - target_grad_cv = grad_cv_operator_target_compiled( + target_grad_cv = grad_cv_operator_target_compiled( + target_fluid_state, time=0.) + if use_av == 1: + target_av_smu = compute_smoothness_compiled( + cv=target_cv, dv=target_fluid_state.dv, grad_cv=target_grad_cv) + elif use_av == 2: + target_grad_t = grad_t_operator_target_compiled( target_fluid_state, time=0.) - if use_av == 1: - target_av_smu = compute_smoothness_compiled( - cv=target_cv, dv=target_fluid_state.dv, grad_cv=target_grad_cv) - elif use_av == 2: - target_grad_t = grad_t_operator_target_compiled( - target_fluid_state, time=0.) - target_av_sbeta = compute_smoothness_beta_compiled( - cv=target_cv, dv=target_fluid_state.dv, grad_cv=target_grad_cv) - target_av_skappa = compute_smoothness_kappa_compiled( - cv=target_cv, dv=target_fluid_state.dv, grad_t=target_grad_t) - target_av_smu = compute_smoothness_mu_compiled( - cv=target_cv, dv=target_fluid_state.dv, grad_cv=target_grad_cv) + target_av_sbeta = compute_smoothness_beta_compiled( + cv=target_cv, dv=target_fluid_state.dv, grad_cv=target_grad_cv) + target_av_skappa = compute_smoothness_kappa_compiled( + cv=target_cv, dv=target_fluid_state.dv, grad_t=target_grad_t) + target_av_smu = compute_smoothness_mu_compiled( + cv=target_cv, dv=target_fluid_state.dv, grad_cv=target_grad_cv) - target_av_smu = force_evaluation(actx, target_av_smu) - target_av_sbeta = force_evaluation(actx, target_av_sbeta) - target_av_skappa = force_evaluation(actx, target_av_skappa) + target_av_smu = force_evaluation(actx, target_av_smu) + target_av_sbeta = force_evaluation(actx, target_av_sbeta) + target_av_skappa = force_evaluation(actx, target_av_skappa) - target_fluid_state = create_fluid_state( - cv=target_cv, temperature_seed=temperature_seed, - smoothness_mu=target_av_smu, smoothness_beta=target_av_sbeta, - smoothness_kappa=target_av_skappa) + target_fluid_state = create_fluid_state( + cv=target_cv, temperature_seed=temperature_seed, + smoothness_mu=target_av_smu, smoothness_beta=target_av_sbeta, + smoothness_kappa=target_av_skappa) ################################## # Set up the boundary conditions # ################################## - from mirgecom.gas_model import project_fluid_state - - def get_target_state_on_boundary(btag): - return project_fluid_state( - dcoll, dd_vol_fluid, - dd_vol_fluid.trace(btag).with_discr_tag(quadrature_tag), - target_fluid_state, gas_model - ) - - flow_ref_state = \ - get_target_state_on_boundary("flow") - - flow_ref_state = force_evaluation(actx, flow_ref_state) - - def _target_flow_state_func(**kwargs): - return flow_ref_state - - flow_boundary = PrescribedFluidBoundary( - boundary_state_func=_target_flow_state_func) - - inflow_ref_state = \ - get_target_state_on_boundary("inflow") - - inflow_ref_state = force_evaluation(actx, inflow_ref_state) - - def _target_inflow_state_func(**kwargs): - return inflow_ref_state - - inflow_boundary = PrescribedFluidBoundary( - boundary_state_func=_target_inflow_state_func) - - outflow_ref_state = \ - get_target_state_on_boundary("outflow") - - outflow_ref_state = force_evaluation(actx, outflow_ref_state) - - def _target_outflow_state_func(**kwargs): - return outflow_ref_state - - outflow_boundary = PrescribedFluidBoundary( - boundary_state_func=_target_outflow_state_func) - #outflow_pressure = 2000 - #outflow_boundary = PressureOutflowBoundary(outflow_pressure) - - if noslip: - if adiabatic: - fluid_wall = AdiabaticNoslipWallBoundary() - else: - fluid_wall = IsothermalWallBoundary(temp_wall) - else: - fluid_wall = AdiabaticSlipBoundary() - - wall_farfield = DirichletDiffusionBoundary(temp_wall) - - if use_injection: - fluid_boundaries = { - flow_bnd.domain_tag: flow_boundary, # pylint: disable=no-member - wall_bnd.domain_tag: fluid_wall # pylint: disable=no-member - } - else: - fluid_boundaries = { - inflow_bnd.domain_tag: inflow_boundary, # pylint: disable=no-member - outflow_bnd.domain_tag: outflow_boundary, # pylint: disable=no-member - inj_bnd.domain_tag: fluid_wall, # pylint: disable=no-member - wall_bnd.domain_tag: fluid_wall # pylint: disable=no-member - } - - wall_boundaries = { - wall_ffld_bnd.domain_tag: wall_farfield # pylint: disable=no-member - } + if init_name == "ACTII": + from y3prediction.actii_y3 import get_boundaries + fluid_boundaries, wall_boundaries, target_boundaries = \ + get_boundaries(dcoll, actx, dd_vol_fluid, dd_vol_wall, + use_injection, quadrature_tag, gas_model, + noslip, adiabatic, temp_wall, target_fluid_state) + elif init_name == "Flash1D": + from y3prediction.flash_utils import get_boundaries + fluid_boundaries, wall_boundaries, target_boundaries = \ + get_boundaries(dcoll, actx, dd_vol_fluid, dd_vol_wall, noslip, + adiabatic, periodic, temp_wall, gas_model, quadrature_tag, + target_fluid_state) # finish initializing the smoothness for non-restarts if not restart_filename: @@ -1961,14 +2075,16 @@ def _get_wall_kappa_fiber(mass, temperature): experimental_kappa(temperature) * puma_kappa(mass_loss_frac) / puma_kappa(0)) - return ( - scaled_insert_kappa * wall_insert_mask - + wall_surround_kappa * wall_surround_mask) + retval = scaled_insert_kappa * wall_insert_mask + if wall_surround_mask is not None: + retval = retval + wall_surround_kappa*wall_surround_mask + return retval def _get_wall_kappa_inert(mass, temperature): - return ( - wall_insert_kappa * wall_insert_mask - + wall_surround_kappa * wall_surround_mask) + retval = wall_insert_kappa * wall_insert_mask + if wall_surround_mask is not None: + retval = retval + wall_surround_kappa*wall_surround_mask + return retval def _get_wall_effective_surface_area_fiber(mass): mass_loss_frac = ( @@ -1989,16 +2105,12 @@ def _mass_loss_rate_fiber(mass, ox_mass, temperature, eff_surf_area): # inert if wall_material == 0: wall_model = WallModel( - heat_capacity=( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask), + heat_capacity=wall_cp, thermal_conductivity_func=_get_wall_kappa_inert) # non-porous elif wall_material == 1: wall_model = WallModel( - heat_capacity=( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask), + heat_capacity=wall_cp, thermal_conductivity_func=_get_wall_kappa_fiber, effective_surface_area_func=_get_wall_effective_surface_area_fiber, mass_loss_func=_mass_loss_rate_fiber, @@ -2006,9 +2118,7 @@ def _mass_loss_rate_fiber(mass, ox_mass, temperature, eff_surf_area): # porous elif wall_material == 2: wall_model = WallModel( - heat_capacity=( - wall_insert_cp * wall_insert_mask - + wall_surround_cp * wall_surround_mask), + heat_capacity=wall_cp, thermal_conductivity_func=_get_wall_kappa_fiber, effective_surface_area_func=_get_wall_effective_surface_area_fiber, mass_loss_func=_mass_loss_rate_fiber, @@ -2173,6 +2283,8 @@ def my_write_viz(step, t, t_wall, fluid_state, wv, wall_kappa, cv = fluid_state.cv dv = fluid_state.dv + mu = fluid_state.viscosity + nu = mu / cv.mass # basic viz quantities, things here are difficult (or impossible) to compute # in post-processing @@ -2231,7 +2343,9 @@ def my_write_viz(step, t, t_wall, fluid_state, wv, wall_kappa, fluid_state.array_context, fluid_state.species_diffusivity ) + cell_Sc = nu / d_alpha_max cell_Pe_mass = char_length_fluid*cv.speed/d_alpha_max + # these are useful if our transport properties # are not constant on the mesh # prandtl @@ -2240,7 +2354,9 @@ def my_write_viz(step, t, t_wall, fluid_state, wv, wall_kappa, viz_ext = [("Re", cell_Re), ("Pe_mass", cell_Pe_mass), - ("Pe_heat", cell_Pe_heat)] + ("Pe_heat", cell_Pe_heat), + ("Sc", cell_Sc)] + fluid_viz_fields.extend(viz_ext) viz_ext = [("char_length_fluid", char_length_fluid), ("char_length_fluid_smooth", smoothed_char_length_fluid)] @@ -2297,7 +2413,7 @@ def my_write_viz(step, t, t_wall, fluid_state, wv, wall_kappa, ("grad_v_x", grad_v[0]), ("grad_v_y", grad_v[1])] if dim == 3: - viz_ext.extend(("grad_v_z", grad_v[2])) + viz_ext.extend([("grad_v_z", grad_v[2])]) viz_ext.extend(("grad_Y_"+species_names[i], grad_y[i]) for i in range(nspecies)) @@ -2667,7 +2783,6 @@ def my_pre_step(step, t, dt, state): from warnings import warn warn("Running gc.collect() to work around memory growth issue " "https://github.com/illinois-ceesd/mirgecom/issues/839") - import gc gc.collect() # Filter *first* because this will be most straightfwd to @@ -2685,7 +2800,9 @@ def my_pre_step(step, t, dt, state): smoothness_beta=av_sbeta, smoothness_kappa=av_skappa) wdv = create_wall_dependent_vars_compiled(wv) - cv = fluid_state.cv # reset cv to limited version + cv = fluid_state.cv + # This re-creation of the state resets *tseed* to current temp + state = make_obj_array([cv, fluid_state.temperature, wv]) try: @@ -2694,7 +2811,6 @@ def my_pre_step(step, t, dt, state): # disable non-constant dt timestepping for now # re-enable when we're ready - do_viz = check_step(step=step, interval=nviz) do_restart = check_step(step=step, interval=nrestart) do_health = check_step(step=step, interval=nhealth) @@ -2881,7 +2997,7 @@ def unfiltered_rhs(t, state): operator_states_quad=operator_fluid_states) """ - ns_rhs, wall_energy_rhs = coupled_ns_heat_operator( + fluid_rhs, wall_energy_rhs = coupled_ns_heat_operator( dcoll=dcoll, gas_model=gas_model, fluid_dd=dd_vol_fluid, wall_dd=dd_vol_wall, @@ -2899,8 +3015,11 @@ def unfiltered_rhs(t, state): chem_rhs = actx.zeros_like(cv) if use_combustion: # conditionals evaluated only once at compile time - chem_rhs = \ - eos.get_species_source_terms(cv, temperature=fluid_state.temperature) + fluid_rhs = ( + fluid_rhs + + eos.get_species_source_terms( + cv, temperature=fluid_state.temperature) + ) ignition_rhs = actx.zeros_like(cv) if use_ignition > 0: @@ -2956,9 +3075,12 @@ def unfiltered_rhs(t, state): sponge_rhs = actx.zeros_like(cv) if use_sponge: - sponge_rhs = _sponge_source(cv=cv) + fluid_rhs = ( + fluid_rhs + + _sponge_source(cv=cv) + ) - fluid_rhs = ns_rhs + chem_rhs + sponge_rhs + ignition_rhs + fluid_rhs = fluid_rhs + chem_rhs + sponge_rhs + ignition_rhs # wall mass loss wall_mass_rhs = actx.zeros_like(wv.mass)