From 5a3ce84fc6bbf8c596850022a3dc9c9b23d2bc4d Mon Sep 17 00:00:00 2001 From: "Eric T. Johnson" Date: Fri, 15 Sep 2023 08:28:47 -0400 Subject: [PATCH] Add new setup for 3D flame propagation through a tube (#2551) --- Exec/science/flame_tube/GNUmakefile | 38 ++ Exec/science/flame_tube/Make.package | 2 + Exec/science/flame_tube/README.md | 5 + Exec/science/flame_tube/_prob_params | 70 +++ .../analysis/vis_3d/andes-single.submit | 34 ++ .../analysis/vis_3d/andes-slice.submit | 17 + .../flame_tube/analysis/vis_3d/andes.submit | 26 + .../analysis/vis_3d/slice_vertical.py | 97 +++ .../analysis/vis_3d/vol-xrb-abar.py | 200 ++++++ .../analysis/vis_3d/vol-xrb-enuc.py | 198 ++++++ .../flame_tube/analysis/vis_3d/vol-xrb.py | 277 +++++++++ Exec/science/flame_tube/initial_model.H | 568 ++++++++++++++++++ .../inputs_He/inputs.He.25cm.static.pslope | 157 +++++ Exec/science/flame_tube/problem_initialize.H | 215 +++++++ .../problem_initialize_state_data.H | 91 +++ Exec/science/flame_tube/problem_tagging.H | 55 ++ 16 files changed, 2050 insertions(+) create mode 100644 Exec/science/flame_tube/GNUmakefile create mode 100644 Exec/science/flame_tube/Make.package create mode 100644 Exec/science/flame_tube/README.md create mode 100644 Exec/science/flame_tube/_prob_params create mode 100644 Exec/science/flame_tube/analysis/vis_3d/andes-single.submit create mode 100644 Exec/science/flame_tube/analysis/vis_3d/andes-slice.submit create mode 100644 Exec/science/flame_tube/analysis/vis_3d/andes.submit create mode 100644 Exec/science/flame_tube/analysis/vis_3d/slice_vertical.py create mode 100644 Exec/science/flame_tube/analysis/vis_3d/vol-xrb-abar.py create mode 100644 Exec/science/flame_tube/analysis/vis_3d/vol-xrb-enuc.py create mode 100644 Exec/science/flame_tube/analysis/vis_3d/vol-xrb.py create mode 100644 Exec/science/flame_tube/initial_model.H create mode 100644 Exec/science/flame_tube/inputs_He/inputs.He.25cm.static.pslope create mode 100644 Exec/science/flame_tube/problem_initialize.H create mode 100644 Exec/science/flame_tube/problem_initialize_state_data.H create mode 100644 Exec/science/flame_tube/problem_tagging.H diff --git a/Exec/science/flame_tube/GNUmakefile b/Exec/science/flame_tube/GNUmakefile new file mode 100644 index 0000000000..04b79a26a3 --- /dev/null +++ b/Exec/science/flame_tube/GNUmakefile @@ -0,0 +1,38 @@ +PRECISION = DOUBLE +PROFILE = FALSE + +DEBUG = FALSE + +DIM = 3 + +COMP = gnu + +USE_MPI = TRUE + +USE_GRAV = TRUE +USE_REACT = TRUE + +USE_ROTATION = FALSE +USE_DIFFUSION = TRUE + +# define the location of the CASTRO top directory +CASTRO_HOME := ../../.. + +USE_JACOBIAN_CACHING = TRUE +USE_CXX_MODEL_PARSER = TRUE +NUM_MODELS := 2 + +# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos +EOS_DIR := helmholtz + +# This sets the network directory in $(MICROPHYSICS_HOME)/networks +NETWORK_DIR := aprox13 + +INTEGRATOR_DIR := VODE + +CONDUCTIVITY_DIR := stellar + +Bpack := ./Make.package +Blocs := . + +include $(CASTRO_HOME)/Exec/Make.Castro diff --git a/Exec/science/flame_tube/Make.package b/Exec/science/flame_tube/Make.package new file mode 100644 index 0000000000..e5cc052427 --- /dev/null +++ b/Exec/science/flame_tube/Make.package @@ -0,0 +1,2 @@ +CEXE_headers += initial_model.H + diff --git a/Exec/science/flame_tube/README.md b/Exec/science/flame_tube/README.md new file mode 100644 index 0000000000..c57c8ba7d5 --- /dev/null +++ b/Exec/science/flame_tube/README.md @@ -0,0 +1,5 @@ +# flame_tube + +This is a slightly modified version of flame_wave that models a 3-d flame +propagating through a tube shaped domain, with a periodic boundary along the +transverse horizontal direction (y). diff --git a/Exec/science/flame_tube/_prob_params b/Exec/science/flame_tube/_prob_params new file mode 100644 index 0000000000..a364fe156e --- /dev/null +++ b/Exec/science/flame_tube/_prob_params @@ -0,0 +1,70 @@ + +dtemp real 3.81e8_rt y + +x_half_max real 1.2e5_rt y + +x_half_width real 3.6e4_rt y + +# cutoff mass fraction of the first species for refinement +X_min real 1.e-4_rt y + +# do we dynamically refine based on density? or based on height? +tag_by_density integer 1 y + +# used for tagging if tag_by_density = 1 +cutoff_density real 500.e0_rt y + +# used if we are refining based on height rather than density +refine_height real 3600 y + +dx_model real 10.0_rt y + +T_hi real 5.e8_rt y + +T_star real 1.e8_rt y + +T_lo real 5.e7_rt y + +dens_base real 2.e6_rt y + +H_star real 500.e0_rt y + +atm_delta real 25.e0_rt y + +fuel1_name character "helium-4" y + +fuel2_name character "" y + +fuel3_name character "" y + +fuel4_name character "" y + +ash1_name character "iron-56" y + +ash2_name character "" y + +ash3_name character "" y + +fuel1_frac real 1.0_rt y + +fuel2_frac real 0.0_rt y + +fuel3_frac real 0.0_rt y + +fuel4_frac real 0.0_rt y + +ash1_frac real 1.0_rt y + +ash2_frac real 0.0_rt y + +ash3_frac real 0.0_rt y + +low_density_cutoff real 1.e-4_rt y + +smallx real 1.e-10_rt y + +x_refine_distance real 1.e30_rt y + +max_hse_tagging_level integer 2 y + +max_base_tagging_level integer 1 y diff --git a/Exec/science/flame_tube/analysis/vis_3d/andes-single.submit b/Exec/science/flame_tube/analysis/vis_3d/andes-single.submit new file mode 100644 index 0000000000..b4dfd409dd --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/andes-single.submit @@ -0,0 +1,34 @@ +#!/bin/bash +#SBATCH -A ast106 +#SBATCH -J flame_tube_vr +#SBATCH -N 1 +#SBATCH -t 4:00:00 + +set -u + +cd $SLURM_SUBMIT_DIR + +source "/ccs/proj/ast106/$USER/mambaforge_$(uname -m)/etc/profile.d/conda.sh" +conda activate andes_yt_dev + +plotfiles=(run_256/flame_tube_25cm_smallplt*) +plotfiles=(run_256_extra_plotfiles/flame_tube_25cm_smallplt0039200) +#plotfile=flame_wave_1000Hz_25cm_smallplt207887 +#plotfile=flame_wave_1000Hz_25cm_smallplt40842 + +#for i in flame_wave_1000Hz_25cm_smallplt*[0-9] +todo=() +for f in "${plotfiles[@]}"; do + dest=${f/run_/analysis_} + # check the last image generated for each plotfile + #if ! [[ -f "${dest}_enuc_annotated_top.png" ]]; then + todo+=("$f") + #fi +done + +export OMP_NUM_THREADS=16 +if [[ ${#todo[@]} -gt 0 ]]; then + srun --cpu-bind=no python ~/dev/Castro/Exec/science/flame_tube/analysis/vis_3d/vol-xrb.py "${todo[@]}" +fi + +echo "done!" diff --git a/Exec/science/flame_tube/analysis/vis_3d/andes-slice.submit b/Exec/science/flame_tube/analysis/vis_3d/andes-slice.submit new file mode 100644 index 0000000000..b26d0d6277 --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/andes-slice.submit @@ -0,0 +1,17 @@ +#!/bin/bash +#SBATCH -A ast106 +#SBATCH -J plots +#SBATCH -N 1 +#SBATCH -t 8:00:00 +#SBATCH -p gpu + +cd $SLURM_SUBMIT_DIR + +#plotfile=flame_wave_1000Hz_25cm_smallplt207887 +#plotfile=flame_wave_1000Hz_25cm_smallplt40842 + +source activate andes_env + +srun python slice_vertical.py flame_wave_1000Hz_25cm_smallplt174562 + + diff --git a/Exec/science/flame_tube/analysis/vis_3d/andes.submit b/Exec/science/flame_tube/analysis/vis_3d/andes.submit new file mode 100644 index 0000000000..6966844303 --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/andes.submit @@ -0,0 +1,26 @@ +#!/bin/bash +#SBATCH -A ast106 +#SBATCH -J plots +#SBATCH -N 1 +#SBATCH -t 8:00:00 + +cd $SLURM_SUBMIT_DIR + +source "/ccs/proj/ast106/$USER/mambaforge_$(uname -m)/etc/profile.d/conda.sh" +conda activate andes_yt_dev + +plotfile=run_256_extra_plotfiles/flame_tube_25cm_smallplt0039200 +#plotfile=flame_wave_1000Hz_25cm_smallplt207887 +#plotfile=flame_wave_1000Hz_25cm_smallplt40842 + +#for i in flame_wave_1000Hz_25cm_smallplt*[0-9] +export OMP_NUM_THREADS=16 +for i in $plotfile; do + dest=${i/run_/analysis_} + if ! [[ -f "${dest}_abar_annotated_top.png" ]]; then + srun --cpu-bind=no python ~/dev/Castro/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-abar.py "${i}" + srun --cpu-bind=no python ~/dev/Castro/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-enuc.py "${i}" + fi +done + +echo "done!" diff --git a/Exec/science/flame_tube/analysis/vis_3d/slice_vertical.py b/Exec/science/flame_tube/analysis/vis_3d/slice_vertical.py new file mode 100644 index 0000000000..10c7ed62e8 --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/slice_vertical.py @@ -0,0 +1,97 @@ +import argparse +import os +import sys + +import matplotlib +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.axes_grid1 import ImageGrid + +import yt +from yt.frontends.boxlib.api import CastroDataset +# assume that our data is in CGS +from yt.units import amu, cm + +matplotlib.use('agg') + + + + +def doit(plotfiles): + + plotfile = sys.argv[1] + ds = CastroDataset(plotfile) + + xmin = ds.domain_left_edge[0] + xmax = ds.domain_right_edge[0] + xctr = 0.5 * (xmin + xmax) + L_x = (2./3.) * (xmax - xmin) + + ymin = ds.domain_left_edge[1] + ymax = ds.domain_right_edge[1] + yctr = 0.5*(ymin + ymax) + + zmin = 0.0*cm + zmax = 1.e4*cm + + zctr = 0.5*(zmin + zmax) + L_z = zmax - zmin + + fig = plt.figure() + + grid = ImageGrid(fig, 111, nrows_ncols=(len(plotfiles), 1), + axes_pad=0.25, label_mode="L", + cbar_mode="single", cbar_size="0.5%") + + for i, pf in enumerate(plotfiles): + + ds = CastroDataset(pf) + + f = "abar" + + sp = yt.SlicePlot(ds, "y", f, origin="native", center=[xctr, yctr, zctr], + width=[L_z, L_x], fontsize="9") + sp.set_buff_size((4800,4800)) + sp.swap_axes() + + sp.set_zlim(f, 4, 5) + sp.set_log(f, False) + sp.set_cmap(f, "plasma_r") + + sp.set_axes_unit("cm") + + sp.annotate_text((0.8, 0.9), "{:5.2f} ms".format(1000.0*float(ds.current_time.in_cgs())), + coord_system="axis", text_args={"color": "black", "size": 9}) + + plot = sp.plots["abar"] + plot.figure = fig + plot.axes = grid[i].axes + plot.cax = grid.cbar_axes[i] + if i < len(plotfiles)-1: + grid[i].axes.xaxis.offsetText.set_visible(False) + + sp._setup_plots() + + fig.set_size_inches(10.0, 3.5) + fig.subplots_adjust(left=0.05, right=0.95, top=0.95, bottom=0.05) + fig.savefig("time_series_3D.pdf") + +if __name__ == "__main__": + + p = argparse.ArgumentParser() + + p.add_argument("--skip", type=int, default=1, + help="interval between plotfiles") + p.add_argument("plotfiles", type=str, nargs="+", + help="list of plotfiles to plot") + + args = p.parse_args() + + plot_prefix = args.plotfiles[0].split("plt")[0] + "plt" + plot_nums = sorted([p.split("plt")[1] for p in args.plotfiles], key=int) + + plotfiles = [] + for n in range(0, len(plot_nums), args.skip): + plotfiles.append("{}{}".format(plot_prefix, plot_nums[n])) + + doit(plotfiles) diff --git a/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-abar.py b/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-abar.py new file mode 100644 index 0000000000..2c87a1349e --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-abar.py @@ -0,0 +1,200 @@ +#!/usr/bin/env python +# pylint: disable=wrong-import-position, wrong-import-order, protected-access +# pylint: disable=too-many-locals, too-many-statements, use-dict-literal + +import matplotlib + +matplotlib.use('agg') + +import sys +from pathlib import Path + +import numpy as np +import yt +from yt.units import cm +from yt.visualization.volume_rendering.api import (PointSource, Scene, + create_volume_source) + +resolution = (1920, 1080) + +# this is for the wdconvect problem + + +def doit(plotfile): + + ds = yt.load(plotfile) + + plotfile = plotfile.replace("run_", "analysis_") + Path(plotfile).parent.mkdir(parents=True, exist_ok=True) + ds._periodicity = (True, True, True) + + field = ('boxlib', 'abar') + ds._get_field_info(field).take_log = False + + sc = Scene() + + + # add a volume: select a sphere + #center = (0, 0, 0) + #R = (5.e8, 'cm') + + #dd = ds.sphere(center, R) + + center = ds.domain_center + # set the center in the vertical direction to be the height of the underlying base layer + center[-1] = 2000*cm + + vol = create_volume_source(ds.all_data(), field=field) + sc.add_source(vol) + + axes_point = min(ds.domain_center) + ref_points = PointSource( + positions=np.array([ + ds.domain_left_edge, + [axes_point, ds.domain_left_edge[1], ds.domain_left_edge[2]], + [ds.domain_left_edge[0], axes_point, ds.domain_left_edge[2]], + [ds.domain_left_edge[0], ds.domain_left_edge[1], axes_point], + # center, + ]), + colors=np.array([ + [1, 1, 1, 0.01], + [1, 0, 0, 0.01], + [0, 1, 0, 0.01], + [0, 0, 1, 0.01], + # [0, 1, 1, 0.01], + ]), + radii=np.array([ + 5, + 3, + 3, + 3, + # 5, + ]) + ) + + # transfer function + vals = [4.25, 4.5, 4.75, 5, 5.25, 5.5] + sigma = 0.05 + + tf = yt.ColorTransferFunction((min(vals), max(vals))) + + tf.clear() + + cmap = "plasma_r" + + for v in vals: + if v < 5.0: + alpha = 0.25 + else: + alpha = 0.75 + + tf.sample_colormap(v, sigma**2, alpha=alpha, colormap=cmap) + + vol.transfer_function = tf + + cam = sc.add_camera(ds, lens_type="perspective") + cam.resolution = resolution + + # view 1 (side) + + setup_side_camera(cam, ds, center) + sc.camera = cam + print(f"side view: {cam!r}", flush=True) + + sc.save(f"{plotfile}_abar_noaxes_side.png", sigma_clip=3.0) + + # set the alpha value for the annotations + alpha = 0.01 + ref_points.colors[:, 3] = alpha + sc.add_source(ref_points) + # sc.annotate_axes(alpha=alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, alpha]) + + sc.save(f"{plotfile}_abar_side.png", sigma_clip=3.0) + + sc.save_annotated(f"{plotfile}_abar_annotated_side.png", + sigma_clip=3.0, label_fmt="%.2f", + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + # view 2 (top) + + # remove the annotation sources for now + print(sc.sources, flush=True) + for key, source in list(sc.sources.items()): + if source is not vol: + sc.sources.pop(key) + + setup_top_camera(cam, ds, center) + sc.camera = cam + print(f"top view: {cam!r}", flush=True) + + sc.save(f"{plotfile}_abar_noaxes_top.png", sigma_clip=3.0) + + # set the alpha value for the annotations + alpha = 0.005 + ref_points.colors[:, 3] = alpha + sc.add_source(ref_points) + # sc.annotate_axes(alpha=alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, alpha]) + + sc.save(f"{plotfile}_abar_top.png", sigma_clip=3.0) + + sc.save_annotated(f"{plotfile}_abar_annotated_top.png", + sigma_clip=3.0, label_fmt="%.2f", + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + +def setup_side_camera(cam, ds, center): + cam.set_position([center[0], + ds.domain_left_edge[1] - ds.domain_width[0] / 2, + center[2]]) # + 0.25 * (ds.domain_right_edge[2] - ds.domain_left_edge[2])] + + # look toward the center + normal = center - cam.position + normal /= np.sqrt(normal.dot(normal)) + + cam.set_focus(center) + cam.switch_orientation(normal_vector=normal, north_vector=[0., 0., 1.]) + # width[0] and width[1] are the length and height of the image plane + # width[2] is the distance from the camera to the image plane + # x is horizontal, z is vertical + # put the image plane at the near side of the domain + width = [ds.domain_width[0], ds.domain_width[2], + min(abs(cam.position[1] - ds.domain_left_edge[1]), + abs(cam.position[1] - ds.domain_right_edge[1]))] + cam.set_width(width) + cam.zoom(1 / 1.1) + + +def setup_top_camera(cam, ds, center): + cam.set_position([ds.domain_center[0], + ds.domain_center[1], + ds.domain_right_edge[2] + ds.domain_width[0] / 2], + north_vector=[0., 1., 0.]) + + normal = center - cam.position + normal /= np.sqrt(normal.dot(normal)) + + cam.switch_orientation(normal_vector=normal, north_vector=[0., 1., 0.]) + # width[0] and width[1] are the length and height of the image plane + # width[2] is the distance from the camera to the image plane + # x is horizontal, y is vertical + # put the image plane at the top of the domain + width = [ds.domain_width[0], ds.domain_width[1], + min(abs(cam.position[2] - ds.domain_left_edge[2]), + abs(cam.position[2] - ds.domain_right_edge[2]))] + cam.set_width(width) + cam.zoom(1 / 1.1) + + +if __name__ == "__main__": + + if not sys.argv[1:]: + sys.exit("ERROR: no plotfile specified") + + for file in sys.argv[1:]: + doit(file) diff --git a/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-enuc.py b/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-enuc.py new file mode 100644 index 0000000000..02bda7f83c --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/vol-xrb-enuc.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python +# pylint: disable=wrong-import-position, wrong-import-order, protected-access +# pylint: disable=too-many-locals, too-many-statements, use-dict-literal + +import matplotlib + +matplotlib.use('agg') + +import sys +from pathlib import Path + +import numpy as np +import yt +from yt.units import cm +from yt.visualization.volume_rendering.api import (PointSource, Scene, + create_volume_source) + +resolution = (1920, 1080) + +# this is for the wdconvect problem + + +def doit(plotfile): + + ds = yt.load(plotfile) + + plotfile = plotfile.replace("run_", "analysis_") + Path(plotfile).parent.mkdir(parents=True, exist_ok=True) + ds._periodicity = (True, True, True) + + field = ('boxlib', 'enuc') + ds._get_field_info(field).take_log = True + + sc = Scene() + + + # add a volume: select a sphere + #center = (0, 0, 0) + #R = (5.e8, 'cm') + + #dd = ds.sphere(center, R) + + center = ds.domain_center + # set the center in the vertical direction to be the height of the underlying base layer + center[-1] = 2000*cm + + vol = create_volume_source(ds.all_data(), field=field) + sc.add_source(vol) + + axes_point = min(ds.domain_center) + ref_points = PointSource( + positions=np.array([ + ds.domain_left_edge, + [axes_point, ds.domain_left_edge[1], ds.domain_left_edge[2]], + [ds.domain_left_edge[0], axes_point, ds.domain_left_edge[2]], + [ds.domain_left_edge[0], ds.domain_left_edge[1], axes_point], + # center, + ]), + colors=np.array([ + [1, 1, 1, 0.01], + [1, 0, 0, 0.01], + [0, 1, 0, 0.01], + [0, 0, 1, 0.01], + # [0, 1, 1, 0.01], + ]), + radii=np.array([ + 5, + 3, + 3, + 3, + # 5, + ]) + ) + + # transfer function + vals = [16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5] + sigma = 0.1 + + tf = yt.ColorTransferFunction((min(vals), max(vals))) + + tf.clear() + + cmap = "viridis" + + for v in vals: + if v < 19.0: + alpha = 0.25 + else: + alpha = 0.75 + + tf.sample_colormap(v, sigma**2, alpha=alpha, colormap=cmap) + + vol.transfer_function = tf + + cam = sc.add_camera(ds, lens_type="perspective") + cam.resolution = resolution + + # view 1 (side) + + setup_side_camera(cam, ds, center) + sc.camera = cam + print(f"side view: {cam!r}", flush=True) + + sc.save(f"{plotfile}_Hnuc_noaxes_side.png") + + # set the alpha value for the annotations + alpha = 0.02 + ref_points.colors[:, 3] = alpha + sc.add_source(ref_points) + # sc.annotate_axes(alpha=alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, alpha]) + + sc.save(f"{plotfile}_Hnuc_side.png") + + sc.save_annotated(f"{plotfile}_Hnuc_annotated_side.png", + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + # view 2 (top) + + # remove the annotation sources for now + print(sc.sources, flush=True) + for key, source in list(sc.sources.items()): + if source is not vol: + sc.sources.pop(key) + + setup_top_camera(cam, ds, center) + sc.camera = cam + print(f"top view: {cam!r}", flush=True) + + sc.save(f"{plotfile}_Hnuc_noaxes_top.png") + + # set the alpha value for the annotations + alpha = 0.005 + ref_points.colors[:, 3] = alpha + sc.add_source(ref_points) + # sc.annotate_axes(alpha=alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, alpha]) + + sc.save(f"{plotfile}_Hnuc_top.png") + + sc.save_annotated(f"{plotfile}_Hnuc_annotated_top.png", + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + +def setup_side_camera(cam, ds, center): + cam.set_position([center[0], + ds.domain_left_edge[1] - ds.domain_width[0] / 2, + center[2]]) # + 0.25 * (ds.domain_right_edge[2] - ds.domain_left_edge[2])] + + # look toward the center + normal = center - cam.position + normal /= np.sqrt(normal.dot(normal)) + + cam.set_focus(center) + cam.switch_orientation(normal_vector=normal, north_vector=[0., 0., 1.]) + # width[0] and width[1] are the length and height of the image plane + # width[2] is the distance from the camera to the image plane + # x is horizontal, z is vertical + # put the image plane at the near side of the domain + width = [ds.domain_width[0], ds.domain_width[2], + min(abs(cam.position[1] - ds.domain_left_edge[1]), + abs(cam.position[1] - ds.domain_right_edge[1]))] + cam.set_width(width) + cam.zoom(1 / 1.1) + + +def setup_top_camera(cam, ds, center): + cam.set_position([ds.domain_center[0], + ds.domain_center[1], + ds.domain_right_edge[2] + ds.domain_width[0] / 2], + north_vector=[0., 1., 0.]) + + normal = center - cam.position + normal /= np.sqrt(normal.dot(normal)) + + cam.switch_orientation(normal_vector=normal, north_vector=[0., 1., 0.]) + # width[0] and width[1] are the length and height of the image plane + # width[2] is the distance from the camera to the image plane + # x is horizontal, y is vertical + # put the image plane at the top of the domain + width = [ds.domain_width[0], ds.domain_width[1], + min(abs(cam.position[2] - ds.domain_left_edge[2]), + abs(cam.position[2] - ds.domain_right_edge[2]))] + cam.set_width(width) + cam.zoom(1 / 1.1) + + +if __name__ == "__main__": + + if not sys.argv[1:]: + sys.exit("ERROR: no plotfile specified") + + for file in sys.argv[1:]: + doit(file) diff --git a/Exec/science/flame_tube/analysis/vis_3d/vol-xrb.py b/Exec/science/flame_tube/analysis/vis_3d/vol-xrb.py new file mode 100644 index 0000000000..f11c403349 --- /dev/null +++ b/Exec/science/flame_tube/analysis/vis_3d/vol-xrb.py @@ -0,0 +1,277 @@ +#!/usr/bin/env python +# pylint: disable=wrong-import-position, wrong-import-order, protected-access, use-dict-literal + +import matplotlib + +matplotlib.use('agg') + +import sys +from pathlib import Path + +import numpy as np +import yt +from yt.units import cm +from yt.visualization.volume_rendering.api import (PointSource, Scene, + create_volume_source) + +resolution = (1920, 1080) + + +def doit(plotfile): + + ds = yt.load(plotfile) + + plotfile = plotfile.replace("run_", "analysis_") + Path(plotfile).parent.mkdir(parents=True, exist_ok=True) + ds._periodicity = (True, True, True) + + # center = ds.domain_center + # # set the center in the vertical direction to be the height of the underlying base layer + # center[-1] = 2000*cm + + axes_point = min(ds.domain_center) + ref_points = PointSource( + positions=np.array([ + ds.domain_left_edge, + [axes_point, ds.domain_left_edge[1], ds.domain_left_edge[2]], + [ds.domain_left_edge[0], axes_point, ds.domain_left_edge[2]], + [ds.domain_left_edge[0], ds.domain_left_edge[1], axes_point], + # center, + ]), + colors=np.array([ + [1, 1, 1, 0.01], + [1, 0, 0, 0.01], + [0, 1, 0, 0.01], + [0, 0, 1, 0.01], + # [0, 1, 1, 0.01], + ]), + radii=np.array([ + 5, + 3, + 3, + 3, + # 5, + ]) + ) + + ds._get_field_info(('boxlib', 'abar')).take_log = False + render_field(ds, ('boxlib', 'abar'), plotfile, + ref_points=ref_points, + side_alpha=0.005, + mid_alpha=0.001, + top_alpha=0.001, + annotated_kwargs={"label_fmt": "%.2f"}, + sigma_clip=3.0) + + ds._get_field_info(('boxlib', 'enuc')).take_log = True + render_field(ds, ('boxlib', 'enuc'), plotfile, + ref_points=ref_points, + side_alpha=0.005, + mid_alpha=0.001, + top_alpha=0.0001, + sigma_clip=3.0) + + del ds + + +def make_tf(field): + if field == ('boxlib', 'abar'): + # transfer function + vals = [4.25, 4.5, 4.75, 5, 5.25, 5.5] + sigma = 0.05 + + tf = yt.ColorTransferFunction((min(vals), max(vals))) + + tf.clear() + + cmap = "plasma_r" + + for v in vals: + if v < 5.0: + alpha = 0.25 + else: + alpha = 0.75 + + tf.sample_colormap(v, sigma**2, alpha=alpha, colormap=cmap) + + elif field == ('boxlib', 'enuc'): + # transfer function + vals = [16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5] + sigma = 0.1 + + tf = yt.ColorTransferFunction((min(vals), max(vals))) + + tf.clear() + + cmap = "viridis" + + for v in vals: + if v < 19.0: + alpha = 0.25 + else: + alpha = 0.75 + + tf.sample_colormap(v, sigma**2, alpha=alpha, colormap=cmap) + + return tf + + +def render_field(ds, field, output_prefix, ref_points, *, + side_alpha=0.01, mid_alpha=0.01, top_alpha=0.01, + annotated_kwargs=None, **kwargs): + + if annotated_kwargs is None: + annotated_kwargs = {} + + output_prefix = f"{output_prefix}_{field[1]}" + + sc = Scene() + + + # add a volume: select a sphere + #center = (0, 0, 0) + #R = (5.e8, 'cm') + + #dd = ds.sphere(center, R) + + vol = create_volume_source(ds.all_data(), field=field) + sc.add_source(vol) + + vol.transfer_function = make_tf(field) + + cam = sc.add_camera(ds, lens_type="perspective") + cam.resolution = resolution + + center = ds.domain_center + # set the center in the vertical direction to be the height of the underlying base layer + center[-1] = 2000*cm + + # view 1 (side) + + setup_side_camera(cam, ds, center) + sc.camera = cam + print(f"side view: {cam!r}", flush=True) + + sc.save(f"{output_prefix}_noaxes_side.png", **kwargs) + + # set the alpha value for the annotations + ref_points.colors[:, 3] = side_alpha + sc.add_source(ref_points) + # sc.annotate_axes(alpha=side_alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, side_alpha]) + + sc.save(f"{output_prefix}_side.png", **kwargs) + + sc.save_annotated(f"{output_prefix}_annotated_side.png", + **kwargs, **annotated_kwargs, + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + remove_annotations(sc, vol) + + # view 2 (upper side) + + setup_mid_camera(cam, ds, center) + sc.camera = cam + print(f"mid view: {cam!r}", flush=True) + + sc.save(f"{output_prefix}_noaxes_mid.png", **kwargs) + + # set the alpha value for the annotations + ref_points.colors[:, 3] = mid_alpha + sc.add_source(ref_points) + # sc.annotate_axes(alpha=mid_alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, mid_alpha]) + + sc.save(f"{output_prefix}_mid.png", **kwargs) + + sc.save_annotated(f"{output_prefix}_annotated_mid.png", + **kwargs, **annotated_kwargs, + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + remove_annotations(sc, vol) + + # view 3 (top) + + setup_top_camera(cam, ds, center) + sc.camera = cam + print(f"top view: {cam!r}", flush=True) + + sc.save(f"{output_prefix}_noaxes_top.png", **kwargs) + + # set the alpha value for the annotations + # ref_points.colors[:, 3] = top_alpha + # sc.add_source(ref_points) + # sc.annotate_axes(alpha=top_alpha, thickness=6) + sc.annotate_domain(ds, color=[1, 1, 1, top_alpha]) + + sc.save(f"{output_prefix}_top.png", **kwargs) + + sc.save_annotated(f"{output_prefix}_annotated_top.png", + **kwargs, **annotated_kwargs, + text_annotate=[[(0.05, 0.05), + f"t = {ds.current_time.d:7.5f} s", + dict(horizontalalignment="left")]]) + + del sc, vol + + +def remove_annotations(sc, vol): + # print(sc.sources, flush=True) + for key, source in list(sc.sources.items()): + if source is not vol: + sc.sources.pop(key) + + +def setup_side_camera(cam, ds, center): + cam.set_position([center[0], + ds.domain_left_edge[1] - ds.domain_width[0], + center[2]], # + 0.25 * (ds.domain_right_edge[2] - ds.domain_left_edge[2])] + north_vector=[0.0, 0.0, 1.0]) + cam.set_focus(center) + + # width[0] and width[1] are the length and height of the image plane + # width[2] is the distance from the camera to the image plane + # x is horizontal, z is vertical + # put the image plane at the near side of the domain + width = [ds.domain_width[0], ds.domain_width[2], + min(abs(cam.position[1] - ds.domain_left_edge[1]), + abs(cam.position[1] - ds.domain_right_edge[1]))] + cam.set_width(width) + cam.zoom(1 / 1.1) + + +def setup_mid_camera(cam, ds, center): + setup_side_camera(cam, ds, center) + + # rotate camera around the x-axis, keeping the distance from the focus constant + cam.rotate(30 * 180.0 / np.pi, rot_center=center, rot_vector=np.array([1.0, 0.0, 0.0])) + + +def setup_top_camera(cam, ds, center): + cam.set_position([ds.domain_center[0], + ds.domain_center[1], + ds.domain_right_edge[2] + ds.domain_width[0]], + north_vector=[0.0, 1.0, 0.0]) + cam.set_focus(center) + + # width[0] and width[1] are the length and height of the image plane + # width[2] is the distance from the camera to the image plane + # x is horizontal, y is vertical + # put the image plane at the top of the domain + width = [ds.domain_width[0], ds.domain_width[1], + min(abs(cam.position[2] - ds.domain_left_edge[2]), + abs(cam.position[2] - ds.domain_right_edge[2]))] + cam.set_width(width) + cam.zoom(1 / 1.1) + + +if __name__ == "__main__": + if not sys.argv[1:]: + sys.exit("ERROR: no plotfiles specified") + + for file in sys.argv[1:]: + doit(file) diff --git a/Exec/science/flame_tube/initial_model.H b/Exec/science/flame_tube/initial_model.H new file mode 100644 index 0000000000..109da68041 --- /dev/null +++ b/Exec/science/flame_tube/initial_model.H @@ -0,0 +1,568 @@ +#ifndef INITIAL_MODEL_H +#define INITIAL_MODEL_H + +#include +#include +#include + +// Create a 1-d hydrostatic, atmosphere with an isothermal region +// (T_star) representing the NS, a hyperbolic tangent rise to a +// peak temperature (T_hi) representing the base of an accreted +// layer, an isoentropic profile down to a lower temperature (T_lo), +// and then isothermal. This can serve as an initial model for a +// nova or XRB. +// +// The temperature profile is: +// +// ^ +// T_hi + ^ . +// | / \ . +// | / \ . +// | / . \ . +// T_star +---------+ \ . +// | . . \ . +// | \ . +// | . . \ . +// T_lo + +----------- . +// | . . . +// +---------+---+---------------> r . +// | \ / +// | atm_delta +// |< H_star>| +// +// +// ^ +// | +// +-- dens_base +// +// dens_base is the density at a height H_star -- just below the rise +// in T up to the peak T_hi. The composition is "ash" in the lower +// isothermal region and "fuel" in the isentropic and upper +// isothermal regions. In the transition region, we apply the same +// hyperbolic tangent profile to interpolate the composition. +// +// The fuel and ash compositions are specified by the fuel?_name, +// fuel?_frac and ash?_name, ash?_frac parameters (name of the species +// and mass fraction). Where ? = 1,2,3. +// +// The model is placed into HSE by the following differencing: +// +// (1/dr) [

_i -

_{i-1} ] = (1/2) [ _i + _{i-1} ] g +// +// This will be iterated over in tandem with the EOS call, +// P(i-1) = P_eos(rho(i-1), T(i-1), X(i-1) +// + +// this version allows for multiple initial models + +using namespace amrex; + +namespace fw { + constexpr Real MAX_ITER = 250; + constexpr Real TOL = 1.e-10_rt; +} + +struct model_t { + + Real xn_base[NumSpec]; + Real xn_star[NumSpec]; + Real xn_perturb[NumSpec]; + + Real dens_base; + Real T_star; + Real T_hi; + Real T_lo; + + Real H_star; + Real atm_delta; + + Real low_density_cutoff; +}; + + +// Evaluate tanh using the exponential form to workaround a PGI bug on Power9 + +AMREX_INLINE +Real evaluate_tanh(const Real z) { + + Real t; + if (std::abs(z) <= 4.0_rt) { + t = (std::exp(z) - std::exp(-z))/(std::exp(z) + std::exp(-z)); + } else if (z < -4.0_rt) { + t = -1.0_rt; + } else { + t = 1.0_rt; + } + + return t; +} + + +AMREX_INLINE +void +generate_initial_model(const int npts_model, const Real xmin, const Real xmax, + const model_t model_params, const int model_num) +{ + + + // Create a 1-d uniform grid that is identical to the mesh that we are + // mapping onto, and then we want to force it into HSE on that mesh. + + // we actually require that the number of points is the same for each + // model, so we'll just set it each time + + model::npts = npts_model; + model::initialized = true; + + if (npts_model > NPTS_MODEL) { + amrex::Error("Error: model has more than NPTS_MODEL points, Increase MAX_NPTS_MODEL"); + } + + // create the grid -- cell centers + + Real dx = (xmax - xmin) / npts_model; + + for (int i = 0; i < npts_model; i++) { + model::profile(model_num).r(i) = + xmin + (static_cast(i) + 0.5_rt) * dx; + } + + + // find the index of the base height + + int index_base = -1; + for (int i = 0; i < npts_model; i++) { + if (model::profile(model_num).r(i) >= xmin + model_params.H_star) { + index_base = i+1; + break; + } + } + + if (index_base == -1) { + amrex::Error("ERROR: invalid base_height"); + } + + + // put the model onto our new uniform grid + + bool fluff = false; + + // determine the conditions at the base -- this is below the atmosphere + + eos_t eos_state; + eos_state.T = model_params.T_star; + eos_state.rho = model_params.dens_base; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model_params.xn_star[n]; + } + + eos(eos_input_rt, eos_state); + + // store the conditions at the base -- we'll use the entropy later + // to constrain the isentropic layer + + Real pres_base = eos_state.p; + + // set an initial temperature profile and composition + + for (int i = 0; i < npts_model; i++) { + + Real xc = model::profile(model_num).r(i) - + (xmin + model_params.H_star) - 1.5_rt * model_params.atm_delta; + + // hyperbolic tangent transition: + + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, model::ispec+n) = + model_params.xn_star[n] + + 0.5_rt * (model_params.xn_base[n] - model_params.xn_star[n]) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // force them to sum to 1 + + Real sumX = 0.0_rt; + for (int n = 0; n < NumSpec; n++) { + sumX += model::profile(model_num).state(i, model::ispec+n); + } + for (int n = 0; n < NumSpec; n++) { + model::profile(model_num).state(i, model::ispec+n) /= sumX; + } + + // temperature profile -- it is constant below the base + + if (i <= index_base) { + model::profile(model_num).state(i, model::itemp) = model_params.T_star; + } else { + model::profile(model_num).state(i, model::itemp) = + model_params.T_star + + 0.5_rt * (model_params.T_hi - model_params.T_star) * + (1.0_rt + evaluate_tanh(xc / (0.5_rt * model_params.atm_delta))); + } + + // the density and pressure will be determined via HSE, + // for now, set them to the base conditions + + model::profile(model_num).state(i, model::idens) = model_params.dens_base; + model::profile(model_num).state(i, model::ipres) = pres_base; + + } + + // make the base thermodynamics consistent for this base point -- that is + // what we will integrate from! + + eos_state.rho = model::profile(model_num).state(index_base, model::idens); + eos_state.T = model::profile(model_num).state(index_base, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(index_base, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + model::profile(model_num).state(index_base, model::ipres) = eos_state.p; + + + // + // HSE + entropy solve + // + + // the HSE state will be done putting creating an isentropic state until + // the temperature goes below T_lo -- then we will do isothermal. + // also, once the density goes below low_density_cutoff, we stop HSE + + bool isentropic = false; + bool flipped = false; // we start out isothermal and then 'flip' to isentropic + + // + // integrate upward + // + + Real entropy_base; + + for (int i = index_base+1; i < npts_model; i++) { + + if ((model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) && !flipped) { + isentropic = true; + flipped = true; + + // now we need to know the entropy we are confining ourselves to + + eos_state.rho = model::profile(model_num).state(i-1, model::idens); + eos_state.T = model::profile(model_num).state(i-1, model::itemp); + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = model::profile(model_num).state(i-1, model::ispec+n); + } + + eos(eos_input_rt, eos_state); + + entropy_base = eos_state.s; + + amrex::Print() << "base density = " << eos_state.rho << " " << eos_state.T << std::endl; + } + + // we've already set initial guesses for density, temperature, and + // composition + + Real dens_zone = model::profile(model_num).state(i, model::idens); + Real temp_zone = model::profile(model_num).state(i, model::itemp); + Real xn[NumSpec]; + for (int n = 0; n < NumSpec; n++) { + xn[n] = model::profile(model_num).state(i, model::ispec+n); + } + + + // + // iteration loop + // + + // start off the Newton loop by saying that the zone has not converged + + bool converged_hse = false; + + Real pres_zone; + Real entropy; + + if (!fluff) { + + Real p_want; + Real drho; + Real dtemp = 0; + + for (int iter = 0; iter < fw::MAX_ITER; iter++) { + + // get the pressure we want from the HSE equation, just the + // zone below the current. Note, we are using an average of + // the density of the two zones as an approximation of the + // interface value -- this means that we need to iterate for + // find the density and pressure that are consistent + + // furthermore, we need to get the entropy that we need, + // which will come from adjusting the temperature in + // addition to the density. + + // HSE differencing + + p_want = model::profile(model_num).state(i-1, model::ipres) + + dx * 0.5_rt * (dens_zone + model::profile(model_num).state(i-1, model::idens)) * gravity::const_grav; + + if (isentropic) { + + // now we have two functions to zero: + // A = p_want - p(rho,T) + // B = entropy_base - s(rho,T) + // We use a two dimensional Taylor expansion and find the deltas + // for both density and temperature + + // now we know the pressure and the entropy that we want, so we + // need to find the temperature and density through a two + // dimensional root find + + // (t, rho) -> (p, s) + + eos_state.T = temp_zone; + eos_state.rho = dens_zone; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = xn[n]; + } + + eos(eos_input_rt, eos_state); + + entropy = eos_state.s; + pres_zone = eos_state.p; + + Real dpt = eos_state.dpdT; + Real dpd = eos_state.dpdr; + Real dst = eos_state.dsdT; + Real dsd = eos_state.dsdr; + + Real A = p_want - pres_zone; + Real B = entropy_base - entropy; + + dtemp = ((dsd / (dpd - 0.5_rt * dx * gravity::const_grav)) * A - B) / + (dsd * dpt / (dpd - 0.5_rt * dx * gravity::const_grav) - dst); + + drho = (A - dpt * dtemp) / (dpd - 0.5_rt * dx * gravity::const_grav); + + dens_zone = amrex::max(0.9_rt * dens_zone, + amrex::min(dens_zone + drho, 1.1_rt * dens_zone)); + + temp_zone = amrex::max(0.9_rt * temp_zone, + amrex::min(temp_zone + dtemp, 1.1_rt * temp_zone)); + + // check if the density falls below our minimum cut-off -- + // if so, floor it + + if (dens_zone < model_params.low_density_cutoff) { + + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } + + // if (A < TOL .and. B < ETOL) then + if (std::abs(drho) < fw::TOL * dens_zone && std::abs(dtemp) < fw::TOL * temp_zone) { + converged_hse = true; + break; + } + + } else { + + // do isothermal + + if (model::profile(model_num).r(i) > xmin + model_params.H_star + 3.0_rt * model_params.atm_delta) { + temp_zone = model_params.T_lo; + } + + // (t, rho) -> (p) + + eos_state.T = temp_zone; + eos_state.rho = dens_zone; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = xn[n]; + } + + eos(eos_input_rt, eos_state); + + entropy = eos_state.s; + pres_zone = eos_state.p; + + Real dpd = eos_state.dpdr; + + drho = (p_want - pres_zone) / (dpd - 0.5_rt * dx * gravity::const_grav); + + dens_zone = amrex::max(0.9_rt * dens_zone, + amrex::min(dens_zone + drho, 1.1_rt * dens_zone)); + + if (std::abs(drho) < fw::TOL * dens_zone) { + converged_hse = true; + break; + } + + if (dens_zone < model_params.low_density_cutoff) { + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + converged_hse = true; + fluff = true; + break; + } + + } + + if (temp_zone < model_params.T_lo) { + temp_zone = model_params.T_lo; + isentropic = false; + } + + } + + if (!converged_hse) { + std::cout << "Error zone " << i << " did not converge in init_1d" << std::endl; + std::cout << "integrate up" << std::endl; + std::cout << dens_zone << " " << temp_zone << std::endl; + std::cout << p_want << " " << entropy_base << " " << entropy << std::endl; + std::cout << drho << " " << dtemp << std::endl; + amrex::Error("Error: HSE non-convergence"); + } + + } else { + // fluff + + dens_zone = model_params.low_density_cutoff; + temp_zone = model_params.T_lo; + } + + + // call the EOS one more time for this zone and then go on to the next + // (t, rho) -> (p) + + eos_state.T = temp_zone; + eos_state.rho = dens_zone; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = xn[n]; + } + + eos(eos_input_rt, eos_state); + + pres_zone = eos_state.p; + + // update the thermodynamics in this zone + + model::profile(model_num).state(i, model::idens) = dens_zone; + model::profile(model_num).state(i, model::itemp) = temp_zone; + model::profile(model_num).state(i, model::ipres) = pres_zone; + + // to make this process converge faster, set the density in the + // next zone to the density in this zone + // model::profile(model_num).state(i+1, model::idens) = dens_zone; + + } + + + // + // integrate down -- using the temperature profile defined above + // + + for (int i = index_base-1; i >= 0; --i) { + + // we already set the temperature and composition profiles + + Real temp_zone = model::profile(model_num).state(i, model::itemp); + Real xn[NumSpec]; + for (int n = 0; n < NumSpec; n++) { + xn[n] = model::profile(model_num).state(i, model::ispec+n); + } + + // use our previous initial guess for density + + Real dens_zone = model::profile(model_num).state(i+1, model::idens); + + // + // iteration loop + // + + // start off the Newton loop by saying that the zone has not converged + + bool converged_hse = false; + + Real pres_zone; + Real p_want; + Real drho; + + for (int iter = 0; iter < fw::MAX_ITER; iter++) { + + // get the pressure we want from the HSE equation, just the + // zone above the current. Note, we are using an average of + // the density of the two zones as an approximation of the + // interface value -- this means that we need to iterate for + // find the density and pressure that are consistent + + // HSE differencing + + p_want = model::profile(model_num).state(i+1, model::ipres) - + dx * 0.5_rt * (dens_zone + model::profile(model_num).state(i+1, model::idens)) * gravity::const_grav; + + // we will take the temperature already defined in gen_model_state + // so we only need to zero: + // A = p_want - p(rho) + + // (t, rho) -> (p) + eos_state.T = temp_zone; + eos_state.rho = dens_zone; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = xn[n]; + } + + eos(eos_input_rt, eos_state); + + pres_zone = eos_state.p; + + Real dpd = eos_state.dpdr; + + Real A = p_want - pres_zone; + + drho = A / (dpd + 0.5_rt * dx * gravity::const_grav); + + dens_zone = amrex::max(0.9_rt * dens_zone, + amrex::min(dens_zone + drho, 1.1_rt * dens_zone)); + + if (std::abs(drho) < fw::TOL * dens_zone) { + converged_hse = true; + break; + } + + } + + if (!converged_hse) { + std::cout << "Error zone " << i << " did not converge in init_1d" << std::endl; + std::cout << "integrate down" << std::endl; + std::cout << dens_zone << " " << temp_zone << std::endl; + std::cout << p_want << std::endl; + std::cout << drho << std::endl; + amrex::Error("Error: HSE non-convergence"); + } + + + // call the EOS one more time for this zone and then go on to the next + // (t, rho) -> (p) + + eos_state.T = temp_zone; + eos_state.rho = dens_zone; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = xn[n]; + } + + eos(eos_input_rt, eos_state); + + pres_zone = eos_state.p; + + // update the thermodynamics in this zone + + model::profile(model_num).state(i, model::idens) = dens_zone; + model::profile(model_num).state(i, model::itemp) = temp_zone; + model::profile(model_num).state(i, model::ipres) = pres_zone; + + } +} +#endif diff --git a/Exec/science/flame_tube/inputs_He/inputs.He.25cm.static.pslope b/Exec/science/flame_tube/inputs_He/inputs.He.25cm.static.pslope new file mode 100644 index 0000000000..24f50beda2 --- /dev/null +++ b/Exec/science/flame_tube/inputs_He/inputs.He.25cm.static.pslope @@ -0,0 +1,157 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 250000000 +stop_time = 3.0 + +# PROBLEM SIZE & GEOMETRY +geometry.coord_sys = 0 # 0 => cart, 1 => RZ 2=>spherical + +geometry.prob_lo = 0 0 0 +geometry.prob_hi = 2.048e5 1.28e4 1.28e4 +amr.n_cell = 1024 64 64 + +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +# 0 = Interior 3 = Symmetry +# 1 = Inflow 4 = SlipWall +# 2 = Outflow 5 = NoSlipWall +# >>>>>>>>>>>>> BC FLAGS <<<<<<<<<<<<<<<< +castro.lo_bc = 2 0 3 +castro.hi_bc = 2 0 2 +geometry.is_periodic = 0 1 0 + +castro.fill_ambient_bc = 1 +castro.ambient_fill_dir = 2 +castro.ambient_outflow_vel = 1 + +castro.domain_is_plane_parallel = 1 + +# WHICH PHYSICS +castro.do_hydro = 1 +castro.do_react = 1 +castro.do_rotation = 0 +castro.do_grav = 1 +castro.do_sponge = 1 + +castro.small_temp = 1.e6 +castro.small_dens = 1.e-5 + +castro.ppm_type = 1 +castro.grav_source_type = 2 +castro.use_pslope = 1 +castro.pslope_cutoff_density = 1.e4 + +gravity.gravity_type = ConstantGrav +gravity.const_grav = -1.5e14 + +#castro.rotational_period = 0.001 +#castro.rotation_include_centrifugal = 0 + +castro.diffuse_temp = 1 +castro.diffuse_cutoff_density_hi = 5.e4 +castro.diffuse_cutoff_density = 2.e4 + +castro.diffuse_cond_scale_fac = 1.0 + +castro.react_rho_min = 1.e2 +castro.react_rho_max = 5.e6 + +castro.react_T_min = 6.e7 + +castro.sponge_upper_density = 1.e2 +castro.sponge_lower_density = 1.e0 +castro.sponge_timescale = 1.e-7 + +# GPU options +castro.hydro_memory_footprint_ratio = 3 + +# TIME STEP CONTROL +castro.cfl = 0.8 # cfl number for hyperbolic system +castro.init_shrink = 0.1 # scale back initial timestep +castro.change_max = 1.1 # max time step growth + +castro.use_retry = 1 +castro.max_subcycles = 16 + +castro.retry_small_density_cutoff = 10.0 + +castro.abundance_failure_tolerance = 0.1 +castro.abundance_failure_rho_cutoff = 1.0 + +castro.state_interp_order = 0 + +# DIAGNOSTICS & VERBOSITY +castro.sum_interval = 10 # timesteps between computing mass +castro.v = 1 # verbosity in Castro.cpp +amr.v = 1 # verbosity in Amr.cpp +amr.run_log_terse = run_log_terse + +# REFINEMENT / REGRIDDING +amr.max_level = 2 # maximum level number allowed +amr.ref_ratio = 4 2 2 2 # refinement ratio +amr.regrid_int = 0 # static grids ftw +amr.blocking_factor = 32 # block factor in grid generation +amr.max_grid_size = 128 +amr.n_error_buf = 2 2 2 2 # number of buffer cells in error est + +# CHECKPOINT FILES +amr.check_file = flame_tube_25cm_chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints +amr.checkpoint_files_output = 1 + +# PLOTFILES +amr.plot_file = flame_tube_25cm_plt # root name of plotfile +amr.plot_per = 2.e-3 # number of seconds between plotfiles +amr.derive_plot_vars = ALL +amr.plot_files_output = 1 + +amr.small_plot_file = flame_tube_25cm_smallplt # root name of plotfile +amr.small_plot_per = 5.e-4 # number of seconds between plotfiles +amr.small_plot_vars = density Temp + +amr.derive_small_plot_vars = abar magvel enuc +amr.file_name_digits = 7 # pad step number with zeros if needed + +# don't write plotfiles when a stop is requested with dump_and_stop +amr.write_plotfile_with_checkpoint = 0 +castro.output_at_completion = 0 + +fab.format = NATIVE_32 + + +# problem initialization + +problem.dtemp = 1.1e9 +# 4.096e4 + prob_lo[0] (1.024e5) +problem.x_half_max = 1.4336e5 +problem.x_half_width = 2048.e0 + +problem.dx_model = 25.e0 + +problem.dens_base = 3.43e6 + +problem.T_star = 3.e8 +problem.T_hi = 3.e8 +problem.T_lo = 8.e6 + +problem.H_star = 2000.e0 +problem.atm_delta = 50.0 + +problem.fuel1_name = "helium-4" +problem.fuel1_frac = 1.0 + +problem.ash1_name = "nickel-56" +problem.ash1_frac = 1.0 + +problem.low_density_cutoff = 1.e-4 + +problem.tag_by_density = 0 +problem.refine_height = 3600.0 +problem.max_base_tagging_level = 3 + +# Microphysics + +integrator.rtol_spec = 1.e-6 +integrator.atol_spec = 1.e-6 + +network.use_tables = 1 + +# vim: ft=conf diff --git a/Exec/science/flame_tube/problem_initialize.H b/Exec/science/flame_tube/problem_initialize.H new file mode 100644 index 0000000000..2b3630584e --- /dev/null +++ b/Exec/science/flame_tube/problem_initialize.H @@ -0,0 +1,215 @@ +#ifndef problem_initialize_H +#define problem_initialize_H + +#include +#include +#include +#include +#include +#include +#include + +AMREX_INLINE +void problem_initialize () +{ + + const Geometry& dgeom = DefaultGeometry(); + + const Real* problo = dgeom.ProbLo(); + const Real* probhi = dgeom.ProbHi(); + + // check to make sure that small_dens is less than low_density_cutoff + // if not, funny things can happen above the atmosphere + + if (castro::small_dens >= 0.99_rt * problem::low_density_cutoff) { + amrex::Error("ERROR: small_dens should be set lower than low_density_cutoff"); + } + + // make sure hse_fixed_temp is the same as T_star, if it's specified + + if (castro::hse_fixed_temp > 0.0_rt && castro::hse_fixed_temp != problem::T_star) { + amrex::Error("ERROR: hse_fixed_temp should be the same as T_star"); + } + + // get the species indices + + bool species_defined = true; + + int ifuel1 = network_spec_index(problem::fuel1_name); + if (ifuel1 < 0) { + species_defined = false; + } + + int ifuel2; + if (!problem::fuel2_name.empty()) { + ifuel2 = network_spec_index(problem::fuel2_name); + if (ifuel2 < 0) { + species_defined = false; + } + } + + int ifuel3; + if (!problem::fuel3_name.empty()) { + ifuel3 = network_spec_index(problem::fuel3_name); + if (ifuel3 < 0) { + species_defined = false; + } + } + + int ifuel4; + if (!problem::fuel4_name.empty()) { + ifuel4 = network_spec_index(problem::fuel4_name); + if (ifuel4 < 0) { + species_defined = false; + } + } + + int iash1 = network_spec_index(problem::ash1_name); + if (iash1 < 0) { + species_defined = false; + } + + int iash2; + if (!problem::ash2_name.empty()) { + iash2 = network_spec_index(problem::ash2_name); + if (iash2 < 0) { + species_defined = false; + } + } + + int iash3; + if (!problem::ash3_name.empty()) { + iash3 = network_spec_index(problem::ash3_name); + if (iash3 < 0) { + species_defined = false; + } + } + + if (! species_defined) { + std::cout << ifuel1 << " " << ifuel2 << " " << ifuel3 << " " << ifuel4 << std::endl; + std::cout << iash1 << " " << iash2 << " "<< iash3 << std::endl; + amrex::Error("ERROR: species not defined"); + } + + model_t model_params; + + // set the composition of the underlying star + + + for (Real &X : model_params.xn_star) { + X = problem::smallx; + } + model_params.xn_star[iash1] = problem::ash1_frac; + if (!problem::ash2_name.empty()) { + model_params.xn_star[iash2] = problem::ash2_frac; + } + if (!problem::ash3_name.empty()) { + model_params.xn_star[iash3] = problem::ash3_frac; + } + + // and the composition of the accreted layer + + for (Real &X : model_params.xn_base) { + X = problem::smallx; + } + model_params.xn_base[ifuel1] = problem::fuel1_frac; + if (!problem::fuel2_name.empty()) { + model_params.xn_base[ifuel2] = problem::fuel2_frac; + } + if (!problem::fuel3_name.empty()) { + model_params.xn_base[ifuel3] = problem::fuel3_frac; + } + if (!problem::fuel4_name.empty()) { + model_params.xn_base[ifuel4] = problem::fuel4_frac; + } + + // check if they sum to 1 + + Real sumX = 0.0_rt; + for (Real X : model_params.xn_star) { + sumX += X; + } + if (std::abs(sumX - 1.0_rt) > NumSpec * problem::smallx) { + amrex::Error("ERROR: ash mass fractions don't sum to 1"); + } + + sumX = 0.0_rt; + for (Real X : model_params.xn_base) { + sumX += X; + } + if (std::abs(sumX - 1.0_rt) > NumSpec * problem::smallx) { + amrex::Error("ERROR: fuel mass fractions don't sum to 1"); + } + + // we are going to generate an initial model from problo(2) to + // probhi(2) with nx_model zones. But to allow for a interpolated + // lower boundary, we'll add 4 ghostcells to this, so we need to + // compute dx + + int nx_model = static_cast((probhi[AMREX_SPACEDIM-1] - problo[AMREX_SPACEDIM-1]) / + problem::dx_model); + + int ng = 4; + + // now generate the initial models + + model_params.dens_base = problem::dens_base; + model_params.T_star = problem::T_star; + model_params.T_hi = problem::T_hi; + model_params.T_lo = problem::T_lo; + + model_params.H_star = problem::H_star; + model_params.atm_delta = problem::atm_delta; + + model_params.low_density_cutoff = problem::low_density_cutoff; + + generate_initial_model(nx_model + ng, + problo[AMREX_SPACEDIM-1] - ng * problem::dx_model, + probhi[AMREX_SPACEDIM-1], + model_params, 0); + + // now create a perturbed model -- we want the same base conditions + // a hotter temperature + + model_params.T_hi = model_params.T_hi + problem::dtemp; + + generate_initial_model(nx_model + ng, + problo[AMREX_SPACEDIM-1] - ng * problem::dx_model, + probhi[AMREX_SPACEDIM-1], + model_params, 1); + + // put the x- and z-center at zero and the y-center in the middle of the domain + + problem::center[0] = 0.0_rt; + problem::center[1] = 0.5_rt * (problo[1] + probhi[1]); + problem::center[2] = 0.0_rt; + + // set the ambient state for the upper boundary condition + + ambient::ambient_state[URHO] = model::profile(0).state(model::npts-1, model::idens); + ambient::ambient_state[UTEMP] = model::profile(0).state(model::npts-1, model::itemp); + for (int n = 0; n < NumSpec; n++) { + ambient::ambient_state[UFS+n] = + ambient::ambient_state[URHO] * model::profile(0).state(model::npts-1, model::ispec+n); + } + + ambient::ambient_state[UMX] = 0.0_rt; + ambient::ambient_state[UMY] = 0.0_rt; + ambient::ambient_state[UMZ] = 0.0_rt; + + // make the ambient state thermodynamically consistent + + eos_t eos_state; + eos_state.rho = ambient::ambient_state[URHO]; + eos_state.T = ambient::ambient_state[UTEMP]; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = ambient::ambient_state[UFS+n] / eos_state.rho; + } + + eos(eos_input_rt, eos_state); + + ambient::ambient_state[UEINT] = eos_state.rho * eos_state.e; + ambient::ambient_state[UEDEN] = eos_state.rho * eos_state.e; +} + +#endif diff --git a/Exec/science/flame_tube/problem_initialize_state_data.H b/Exec/science/flame_tube/problem_initialize_state_data.H new file mode 100644 index 0000000000..30a9ad8e64 --- /dev/null +++ b/Exec/science/flame_tube/problem_initialize_state_data.H @@ -0,0 +1,91 @@ +#ifndef problem_initialize_state_data_H +#define problem_initialize_state_data_H + +#include +#include +#include +#include +#include +#include +#include + +using namespace amrex; + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_initialize_state_data (int i, int j, int k, + Array4 const& state, + const GeometryData& geomdata) +{ + + const Real* dx = geomdata.CellSize(); + const Real* problo = geomdata.ProbLo(); + + Real x = problo[0] + dx[0] * (static_cast(i) + 0.5_rt) - problem::center[0]; + + Real z = problo[2] + dx[2] * (static_cast(k) + 0.5_rt); + + // lateral distance and height + + Real r; + Real height; + + // we want the flame front to be parallel to the y-axis + r = x; + height = z; + + // blending factor + + Real f; + + if (r < problem::x_half_max) { + f = 1.0_rt; + + } else if (r > problem::x_half_max + problem::x_half_width) { + f = 0.0_rt; + + } else { + f = -(r - problem::x_half_max) / problem::x_half_width + 1.0_rt; + } + + state(i,j,k,URHO) = f * interpolate(height, model::idens, 1) + + (1.0_rt - f) * interpolate(height, model::idens, 0); + + state(i,j,k,UTEMP) = f * interpolate(height, model::itemp, 1) + + (1.0_rt - f) * interpolate(height, model::itemp, 0); + + Real temppres = f * interpolate(height, model::ipres, 1) + + (1.0_rt - f) * interpolate(height, model::ipres, 0); + + for (int n = 0; n < NumSpec; n++) { + state(i,j,k,UFS+n) = f * interpolate(height, model::ispec+n, 1) + + (1.0_rt - f) * interpolate(height, model::ispec+n, 0); + } + + eos_t eos_state; + eos_state.rho = state(i,j,k,URHO); + eos_state.T = state(i,j,k,UTEMP); + eos_state.p = temppres; + for (int n = 0; n < NumSpec; n++) { + eos_state.xn[n] = state(i,j,k,UFS+n); + } + + eos(eos_input_rp, eos_state); + + state(i,j,k,UTEMP) = eos_state.T; + state(i,j,k,UEINT) = eos_state.rho * eos_state.e; + state(i,j,k,UEDEN) = state(i,j,k,UEINT); + + // Initial velocities = 0 + + state(i,j,k,UMX) = 0.e0_rt; + state(i,j,k,UMY) = 0.e0_rt; + state(i,j,k,UMZ) = 0.e0_rt; + + // convert to partial densities + + for (int n = 0; n < NumSpec; n++) { + state(i,j,k,UFS+n) = state(i,j,k,URHO) * state(i,j,k,UFS+n); + } +} + +#endif diff --git a/Exec/science/flame_tube/problem_tagging.H b/Exec/science/flame_tube/problem_tagging.H new file mode 100644 index 0000000000..76ef828dcd --- /dev/null +++ b/Exec/science/flame_tube/problem_tagging.H @@ -0,0 +1,55 @@ +#ifndef problem_tagging_H +#define problem_tagging_H +#include +#include +#include + +/// +/// Define problem-specific tagging criteria +/// +/// @param i x-index +/// @param j y-index +/// @param k z-index +/// @param tag tag array (TagBox) +/// @param state simulation state (Fab) +/// @param level AMR level +/// @param geomdata geometry data +/// +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void problem_tagging(int i, int j, int k, + Array4 const& tag, + Array4 const& state, + int level, const GeometryData& geomdata) +{ + + GpuArray loc; + position(i, j, k, geomdata, loc); + + if (problem::tag_by_density) { + if (state(i,j,k,URHO) > problem::cutoff_density && + state(i,j,k,UFS) / state(i,j,k,URHO) > problem::X_min) { + + Real dist = std::abs(loc[0] - problem::center[0]); + if (level < problem::max_hse_tagging_level && dist < problem::x_refine_distance) { + tag(i,j,k) = TagBox::SET; + } + } + + if (state(i,j,k,URHO) > problem::cutoff_density) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + + } else { + + // tag everything below a certain height + if (loc[AMREX_SPACEDIM-1] < problem::refine_height) { + if (level < problem::max_base_tagging_level) { + tag(i,j,k) = TagBox::SET; + } + } + } + +} +#endif