Skip to content

Commit

Permalink
Merge branch 'development' into ppm_one_sided_second
Browse files Browse the repository at this point in the history
  • Loading branch information
zingale authored Nov 17, 2024
2 parents 68702e3 + 2fb7228 commit c8ecaa1
Show file tree
Hide file tree
Showing 28 changed files with 1,067 additions and 246 deletions.
21 changes: 21 additions & 0 deletions CHANGES.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,24 @@
# 24.11

* a new well-balanced method was added to the CTU PPM solver. This
does the characteristic projection only on the perturbed pressure
and then adds back in the hydrostatic pressure. It can be enabled
via `castro.ppm_well_balanced` (#2945)

* fixed a bug in the div{U} calculation for artificial viscosity on
symmetry boundaries (#2983)

* more development on 2D spherical geometry (#2973, #2975, #2981)

* updates to the massive star plotting scripts (#2979)

* add some new checks to prevent running unsupported combinations
of solvers (#2978)

* documentation updates (#2977)

* `flame_wave` can now be run in 1D (#2976)

# 24.10

* update initial model for `subchandra` when doing ASE NSE (#2970)
Expand Down
52 changes: 26 additions & 26 deletions Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out
Original file line number Diff line number Diff line change
@@ -1,40 +1,40 @@
plotfile = plt_00025
time = 0.00026083835470476599
time = 0.00026084061812361401
variables minimum value maximum value
density 5.4596904436e-13 1.2720098871e-12
xmom 1.2327065163e-07 1.4587103381e-07
density 5.4596904436e-13 1.2717694091e-12
xmom 1.2328047524e-07 1.4591226172e-07
ymom 0 0
zmom 0 0
rho_E 0.021940622286 0.039199658391
rho_e 0.0068091599527 0.032333529333
Temp 100.00003116 217.3564598
rho_X 5.4596904436e-13 1.2720098871e-12
rad 7.5662679486e-07 1.4083425328e-05
pressure 0.0045394399687 0.021555686223
kineng 0.0064467244856 0.015131462333
soundspeed 117717.6283 173551.23637
rho_E 0.021940622285 0.039194728059
rho_e 0.0068091599514 0.032333527949
Temp 100.00003115 217.29660609
rho_X 5.4596904436e-13 1.2717694091e-12
rad 7.5662675907e-07 1.4083418723e-05
pressure 0.0045394399678 0.0215556853
kineng 0.0064473434191 0.015131462334
soundspeed 117717.62829 173527.33922
Gamma_1 1.6666666667 1.6666666667
MachNumber 0.60689740202 1.9999997205
uplusc 269483.97988 362284.15407
uminusc -66810.704217 117717.5954
entropy 2458725989.8 2507243428
MachNumber 0.6068973935 1.9999997207
uplusc 269518.17809 362247.33487
uminusc -66804.446805 117717.59541
entropy 2458725989.8 2507184972.2
magvort 0 0
divu -20286.061564 245.50876324
eint_E 12471696011 27108028479
eint_e 12471696011 27108028479
logden -12.26283198 -11.895509513
StateErr_0 5.4596904436e-13 1.2720098871e-12
StateErr_1 100.00003116 217.3564598
divu -20283.836079 245.94831837
eint_E 12471696009 27100563708
eint_e 12471696009 27100563708
logden -12.26283198 -11.895591626
StateErr_0 5.4596904436e-13 1.2717694091e-12
StateErr_1 100.00003115 217.29660609
StateErr_2 1 1
X(X) 1 1
abar 1 1
x_velocity 102258.95554 235435.2237
x_velocity 102299.57948 235435.22371
y_velocity 0 0
z_velocity 0 0
magvel 102258.95554 235435.2237
radvel -235435.2237 114038.34336
circvel 0 0.005524271728
magmom 1.2327065163e-07 1.4587103381e-07
magvel 102299.57948 235435.22371
radvel -235435.22371 114042.93363
circvel 0 0.0047841596539
magmom 1.2328047524e-07 1.4591226172e-07
angular_momentum_x 0 0
angular_momentum_y 0 0
angular_momentum_z -0 -0
Expand Down
40 changes: 40 additions & 0 deletions Exec/science/xrb_spherical/GNUmakefile
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
PRECISION = DOUBLE
PROFILE = FALSE

DEBUG = FALSE

DIM = 2

COMP = gnu

USE_MPI = TRUE

USE_GRAV = TRUE
USE_REACT = FALSE

USE_ROTATION = FALSE
USE_DIFFUSION = FALSE

# define the location of the CASTRO top directory
CASTRO_HOME ?= ../../..

USE_JACOBIAN_CACHING = TRUE
USE_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 := subch_base

INTEGRATOR_DIR := VODE

CONDUCTIVITY_DIR := stellar

PROBLEM_DIR ?= ./

Bpack := $(PROBLEM_DIR)/Make.package
Blocs := $(PROBLEM_DIR)

include $(CASTRO_HOME)/Exec/Make.Castro
2 changes: 2 additions & 0 deletions Exec/science/xrb_spherical/Make.package
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
CEXE_headers += initial_model.H

6 changes: 6 additions & 0 deletions Exec/science/xrb_spherical/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
# xrb_spherical

This is the full-star XRB flame setup based on flame_wave.
This setup uses a spherical 2D geometry to model XRB flame
on a spherical shell with initial temperature perturbation
on the north pole.
68 changes: 68 additions & 0 deletions Exec/science/xrb_spherical/_prob_params
Original file line number Diff line number Diff line change
@@ -0,0 +1,68 @@

dtemp real 3.81e8_rt y

theta_half_max real 1.745e-2_rt y

theta_half_width real 4.9e-3_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

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 string "helium-4" y

fuel2_name string "" y

fuel3_name string "" y

fuel4_name string "" y

ash1_name string "iron-56" y

ash2_name string "" y

ash3_name string "" 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

r_refine_distance real 1.e30_rt y

max_hse_tagging_level integer 2 y

max_base_tagging_level integer 1 y
56 changes: 56 additions & 0 deletions Exec/science/xrb_spherical/analysis/r_profile.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
#!/usr/bin/env python3

# Spherical R profile at different theta

import os
import sys
import yt
import matplotlib.pyplot as plt
import numpy as np
from functools import reduce
import itertools

import matplotlib.ticker as ptick
from yt.frontends.boxlib.api import CastroDataset
from yt.units import cm


plotfile = sys.argv[1]
ds = CastroDataset(plotfile)

rmin = ds.domain_left_edge[0]
rmax = rmin + 5000.0*cm
#rmax = ds.domain_right_edge[0]
print(ds.domain_left_edge[1])
fig, _ax = plt.subplots(2,2)

axes = list(itertools.chain(*_ax))

fig.set_size_inches(7.0, 8.0)

fields = ["Temp", "density", "x_velocity", "y_velocity"]
nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"]

# 4 rays at different theta values
thetal = ds.domain_left_edge[1]
thetar = ds.domain_right_edge[1]
thetas = [thetal, 0.25*thetar, 0.5*thetar, 0.75*thetar]

for i, f in enumerate(fields):

for theta in thetas:
# simply go from (rmin, theta) -> (rmax, theta). Doesn't need to convert to physical R-Z
ray = ds.ray((rmin, theta, 0*cm), (rmax, theta, 0*cm))
isrt = np.argsort(ray["t"])
axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(float(theta)))

axes[i].set_xlabel(r"$r$ (cm)")
axes[i].set_ylabel(nice_names[i])
axes[i].set_yscale("symlog")

if i == 0:
axes[0].legend(frameon=False, loc="lower left")

#fig.set_size_inches(10.0, 9.0)
plt.tight_layout()
plt.savefig("{}_profiles.png".format(os.path.basename(plotfile)))
94 changes: 94 additions & 0 deletions Exec/science/xrb_spherical/analysis/slice.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,94 @@
#!/usr/bin/env python3

import sys
import os
import yt
import numpy as np
import matplotlib.pyplot as plt
from yt.frontends.boxlib.api import CastroDataset

from yt.units import cm

"""
Given a plot file and field name, it gives slice plots at the top,
middle, and bottom of the domain (shell).
"""

def slice(fname:str, field:str,
loc: str = "top", width_factor: float = 3.0) -> None:
"""
A slice plot of the dataset for Spherical2D geometry.
Parameter
=======================
fname: plot file name
field: field parameter
loc: location on the domain. {top, mid, bot}
"""

ds = CastroDataset(fname)
currentTime = ds.current_time.in_units("s")
print(f"Current time of this plot file is {currentTime} s")

# Some geometry properties
rr = ds.domain_right_edge[0].in_units("km")
rl = ds.domain_left_edge[0].in_units("km")
dr = rr - rl
r_center = 0.5 * (rr + rl)

thetar = ds.domain_right_edge[1]
thetal = ds.domain_left_edge[1]
dtheta = thetar - thetal
theta_center = 0.5 * (thetar + thetal)

# Domain width of the slice plot
width = width_factor * dr
box_widths = (width, width)

loc = loc.lower()
loc_options = ["top", "mid", "bot"]

if loc not in loc_options:
raise Exception("loc parameter must be top, mid or bot")

# Centers for the Top, Mid and Bot panels
centers = {"top":(r_center*np.sin(thetal)+0.5*width, r_center*np.cos(thetal)),
"mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)),
"bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))}

# Note we can also set center during SlicePlot, however then we would enter in [r_center, theta_center, 0]
# rather than the physical R-Z coordinate if we do it via sp.set_center

sp = yt.SlicePlot(ds, 'phi', field, width=box_widths)
sp.set_center(centers[loc])

sp.set_cmap(field, "viridis")
if field in ["x_velocity", "y_velocity", "z_velocity"]:
sp.set_cmap(field, "coolwarm")
elif field == "Temp":
sp.set_zlim(f, 5.e7, 2.5e9)
sp.set_cmap(f, "magma_r")
elif field == "enuc":
sp.set_zlim(f, 1.e18, 1.e20)
elif field == "density":
sp.set_zlim(f, 1.e-3, 5.e8)

# sp.annotate_text((0.05, 0.05), f"{currentTime.in_cgs():8.5f} s")
sp.save(f"{ds}_{loc}")

if __name__ == "__main__":

if len(sys.argv) < 3:
raise Exception("Please enter parameters in order of: fname field_name width_factor[optional] loc[optional]")

fname = sys.argv[1]
field = sys.argv[2]
loc = "top"
width_factor = 3.0

if len(sys.argv) == 4:
width_factor = float(sys.argv[3])
elif len(sys.argv) > 4:
loc = sys.argv[4]

slice(fname, field, loc=loc, width_factor=width_factor)
1 change: 1 addition & 0 deletions Exec/science/xrb_spherical/initial_model.H
Loading

0 comments on commit c8ecaa1

Please sign in to comment.