Skip to content

Commit

Permalink
Added working radiation tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@rbooth200
rbooth200 committed Jul 5, 2023
1 parent ac1aa1b commit a65382b
Show file tree
Hide file tree
Showing 7 changed files with 369 additions and 81 deletions.
85 changes: 85 additions & 0 deletions test_files/irradiation1.par
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
#
# Parameterfile for atmosphere simulations
# Do not use tabs, only white spaces as value separators
#
SPECIES_FILE irradiation.spc
DENSITY_FLOOR 1e-20
PARI_TIME_TMAX 1.01e11
CFL_BREAK_TIME 1e3
PARI_TIME_OUTPUT 2.e10
PARI_CFLFACTOR 0.1
PARI_TIME_DT 1e10
PARI_DPHI_FACTOR 1.5e-2
ENERGY_EPSILON 0.01
MAX_TIMESTEP_CHANGE 1.01
XI_RAD 2.
# Control parameters
#
# Problem numbers explained: shock_tube, no grav == 1, planet with grav == 2#
# Friction solver: analytic == 0, numerical == 1
# Collision model: constant == C, physical == P
PARI_PROBLEM_NUMBER 2
PARI_NUM_SPECIES 1
PARI_NUM_BANDS 2
NUM_BANDS_OUT 1
DO_HYDRO 0
#CONST_TEMP 1430.
PARI_CONST_TEMP 2.7
INIT_TEMPERATURE_MODEL P
INIT_J_FACTOR 1e0
INIT_T_TEMP 2.7
USE_PLANET_TEMPERATURE 1
PARI_TPLANET 350.
PARI_CORE_CV 1e9

# Solar radiation temps at 0.1 AU: 645K for bond albedo=0.75, 913K for bond albedo = 0
# print 5777*(1./20./2.)**0.5*(1-0.75)**0.25

FRICTION_SOLVER 1
PARI_COLL_MODEL P
PARI_ALPHA_COLL 1.
PARI_ORDER 2
PARI_USE_RADIATION 1
PARI_OPACITY_MODEL C
PRESSURE_BROADENING 0.
CONSTOPA_SOLAR_FACTOR 0.35
PARI_CONST_OPAC 1.0

PARI_SELF_GRAV_SWITCH 0
PARI_LINEAR_GRAV 0
PARI_DEBUGLEVEL 1
PARI_INIT_STATIC 1
USE_TIDES 1
#Physical parameters

PARI_PLANET_MASS 224.
PARI_INIT_DATA_SWITCH 1
PARI_INIT_DATA_U1 1e-5
PARI_INIT_DATA_U2 1.0
PARI_INIT_DATA_U3 1e-68

#Radiation parameters
PARI_LAM_MIN 9.0e-2
PARI_LAM_MAX 9.0e-2
PARI_LAM_PER_DECADE 0.1
PARI_TSTAR 6070.
PARI_RSTAR 0.1227
PARI_UVSTAR 1.e1

PARI_PLANET_DIST 0.05

#Numerical parameters
PARI_BOUND_TYPE 2
PARI_BOUND_TYPE_LEFT 2
PARI_BOUND_TYPE_RIGHT 1
PARI_DOMAIN_DX 1e-3
PARI_GRID_TYPE 1
PARI_GEOMETRY 2
PARI_CELLS_PER_DECADE 1000
#Domain min/max are 8 me core radius for earth density, and 100*hill radius
PARI_DOMAIN_MIN 9.200e9
PARI_DOMAIN_MAX 9.2e11

#Scenarios
PARI_INIT_WIND 1
REVERSE_HYDROSTAT_CONSTRUCTION 1
85 changes: 85 additions & 0 deletions test_files/irradiation2.par
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
#
# Parameterfile for atmosphere simulations
# Do not use tabs, only white spaces as value separators
#
SPECIES_FILE mix_static_paperplot.spc
DENSITY_FLOOR 1e-20
PARI_TIME_TMAX 1.01e13
CFL_BREAK_TIME 1e3
PARI_TIME_OUTPUT 2.e13
PARI_CFLFACTOR 0.1
PARI_TIME_DT 1e10
PARI_DPHI_FACTOR 1.5e-2
ENERGY_EPSILON 0.01
MAX_TIMESTEP_CHANGE 1.01
XI_RAD 2.
# Control parameters
#
# Problem numbers explained: shock_tube, no grav == 1, planet with grav == 2#
# Friction solver: analytic == 0, numerical == 1
# Collision model: constant == C, physical == P
PARI_PROBLEM_NUMBER 2
PARI_NUM_SPECIES 1
PARI_NUM_BANDS 2
NUM_BANDS_OUT 1
DO_HYDRO 0
#CONST_TEMP 1430.
PARI_CONST_TEMP 2.7
INIT_TEMPERATURE_MODEL P
INIT_J_FACTOR 1e0
INIT_T_TEMP 2.7
USE_PLANET_TEMPERATURE 1
PARI_TPLANET 350.
PARI_CORE_CV 1e9

# Solar radiation temps at 0.1 AU: 645K for bond albedo=0.75, 913K for bond albedo = 0
# print 5777*(1./20./2.)**0.5*(1-0.75)**0.25

FRICTION_SOLVER 1
PARI_COLL_MODEL P
PARI_ALPHA_COLL 1.
PARI_ORDER 2
PARI_USE_RADIATION 1
PARI_OPACITY_MODEL C
PRESSURE_BROADENING 0.
CONSTOPA_SOLAR_FACTOR 0.35
PARI_CONST_OPAC 1.0

PARI_SELF_GRAV_SWITCH 0
PARI_LINEAR_GRAV 0
PARI_DEBUGLEVEL 1
PARI_INIT_STATIC 1
USE_TIDES 1
#Physical parameters

PARI_PLANET_MASS 224.
PARI_INIT_DATA_SWITCH 1
PARI_INIT_DATA_U1 1e-5
PARI_INIT_DATA_U2 1.0
PARI_INIT_DATA_U3 1e-68

#Radiation parameters
PARI_LAM_MIN 9.0e-2
PARI_LAM_MAX 9.0e-2
PARI_LAM_PER_DECADE 0.1
PARI_TSTAR 6070.
PARI_RSTAR 0.
PARI_UVSTAR 0.e1

PARI_PLANET_DIST 0.05

#Numerical parameters
PARI_BOUND_TYPE 2
PARI_BOUND_TYPE_LEFT 2
PARI_BOUND_TYPE_RIGHT 1
PARI_DOMAIN_DX 1e-3
PARI_GRID_TYPE 1
PARI_GEOMETRY 2
PARI_CELLS_PER_DECADE 1000
#Domain min/max are 8 me core radius for earth density, and 100*hill radius
PARI_DOMAIN_MIN 9.200e9
PARI_DOMAIN_MAX 9.2e11

#Scenarios
PARI_INIT_WIND 1
REVERSE_HYDROSTAT_CONSTRUCTION 1
76 changes: 41 additions & 35 deletions test_files/load_aiolos.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,14 +8,11 @@ def load_aiolos_snap(filename):
('pressure', 'f8'),
('velocity', 'f8'),
('temperature', 'f8'),
('dt_cs', 'f8'),
('dt_v', 'f8'),
('soundspeed', 'f8'),
('dt_E', 'f8'),


]
cols=[0,1,2,3,10,11,12,13,14,15,16]
('soundspeed', 'f8'),
('potential', 'f8'),
('cooling', 'f8'),
('heating', 'f8')]
cols=[0,1,2,3,10,11,12,15,19,21,22]

data = np.genfromtxt(filename, dtype=dtype, usecols=cols)
return data
Expand All @@ -34,14 +31,14 @@ def load_aiolos_params(filename):
params[key.strip()] = val.strip()
return params

def load_aiolos_diag(filename, bands=None, species=None):
if bands is None or species is None:
def load_aiolos_diag(filename, bands_in=None, bands_out=1, species=None):
if bands_in is None or species is None:
# Work out the number of bands:
with open(filename) as f:
line = f.readline()
cols = len(line.strip().split())

def c(b, s): return 5 + 7*b + 3*b*s + s
def c(b, s): return 5 + (3 + s)*b + s + (4+2*s)*bands_out
res = []
b = 1
while True:
Expand All @@ -58,42 +55,45 @@ def c(b, s): return 5 + 7*b + 3*b*s + s
raise ValueError("Error when finding bands and species")
if len(res) > 1:
for b, s in res:
if b == bands or s == species:
bands = b
if b == bands_in or s == species:
bands_in = b
species = s
break
else:
raise ValueError("Error determining the number of bands and"
"species, help out by specifying them")
else:
bands, species = res[0]
bands_in, species = res[0]

# Construct the dtype:
def add_cols(name, dt='f8'):
def add_cols(name, bands=bands_out, dt='f8'):
return [ (name + str(i), dt) for i in range(bands)]
def add_kappa_cols():
cols = []
for i in range(bands):
for i in range(bands_in):
for j in range(species):
cols += [('kappa_P(Tsun)' + str(i) + '_' + str(j), 'f8')]
for i in range(bands_out):
for j in range(species):
cols += [('kappa_R' + str(i) + '_' + str(j), 'f8')]
for j in range(species):
cols += [('kappa_P' + str(i) + '_' + str(j), 'f8')]
cols += [('kappa_P(Tsun)' + str(i) + '_' + str(j), 'f8')]
return cols

dtype = [('x', 'f8'),
('Jtot', 'f8'),
('Stot', 'f8')]
dtype += add_cols('J')
dtype += add_cols('S')
dtype += add_cols('dS')
dtype += add_cols('rho_kappaR')
dtype += add_cols('tau_cell')
dtype += add_cols('tau_radial')
dtype += add_cols('J',bands_out)
dtype += add_cols('S',bands_in)
dtype += add_cols('dS', bands_in)
dtype += add_cols('rho_kappaR',bands_out)
dtype += add_cols('tau_cell',bands_out)
dtype += add_cols('tau_radial', bands_in)
dtype += add_kappa_cols()
dtype += [('T_rad', 'f8')]
dtype += [ ('T_gas' + str(i), 'f8') for i in range(species)]
dtype += [('pressure' , 'f8')]
dtype += add_cols('Flux')
dtype += add_cols('Flux',bands_out)

# Load the data
return np.genfromtxt(filename, dtype=dtype)
Expand Down Expand Up @@ -125,9 +125,11 @@ def load_aiolos_species(filename):
import sys
import matplotlib.pyplot as plt

f, ax = plt.subplots(3,1)
f, ax = plt.subplots(4,1, sharex=True)

for snap in sys.argv[1:]:
snaps = [s for s in sys.argv[1:] if s[0] == 'o'] + [s for s in sys.argv[1:] if s[0] != 'o']

for snap in snaps:
if snap.startswith('output'):
data = load_aiolos_snap(snap)
if 'Gas' in snap:
Expand All @@ -141,26 +143,30 @@ def load_aiolos_species(filename):

ax[1].semilogx(data['x'], data['temperature'], ls='',
marker=m, label=snap)
#ax[1].semilogx(data['x'], data['pressure']/data['density']**1.29, ls='',
# marker=m, label=snap)

#ax[2].semilogx(data['x'], data['dt_cs'], ls='-',
# marker=m, label=snap)
#ax[2].semilogx(data['x'], data['dt_v'], ls='--',
# marker=m, label=snap)
#ax[2].semilogx(data['x'], data['dt_E'], ls=':',
# marker=m, label=snap)
ax[2].semilogx(data['x'], data['velocity'], ls='',
marker=m, label=snap)
ax[2].semilogx(data['x'], data['soundspeed'], ls='--',
marker=m)

l,= ax[3].loglog(data['x'], data['heating'], ls='-',
label=snap)
ax[3].semilogx(data['x'], np.abs(data['cooling']), ls='--',
c=l.get_color())
#ax[2].semilogx(data['x'], data['soundspeed'], ls='--',
# marker=m)
elif snap.startswith('diagnostic'):
data = load_aiolos_diag(snap)
data = load_aiolos_diag(snap, species=2)
ax[1].semilogx(data['x'], data['T_rad'], ls='',
marker='^', label=snap)
#ax[3].semilogx(data['x'], np.abs(data['dS0']), marker='^', ls='')

ax[-1].set_xlabel('r')
ax[0].set_ylabel('density')
ax[1].set_ylabel('temperature')
ax[1].set_ylim(0, 4000)
ax[2].set_ylabel('velocity')
ax[3].set_ylabel('heating / cooling')


ax[1].legend(ncol=3)
Expand Down
56 changes: 28 additions & 28 deletions test_files/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,42 +37,42 @@ bool run_test(std::string test_name, std::string spc = "hydrogen.spc",

int main()
{
int count = 0 ;
count += !run_test("shock_tube1.par") ;
count += !run_test("shock_tube2.par") ;
count += !run_test("shock_tube3.par") ;
count += !run_test("shock_tube4.par") ;
count += !run_test("shock_tube5.par") ;
count += !run_test("shock_tube6.par") ;
count += !run_test("shock_tube7.par") ;
int fail_count = 0 ;
fail_count += !run_test("shock_tube1.par") ;
fail_count += !run_test("shock_tube2.par") ;
fail_count += !run_test("shock_tube3.par") ;
fail_count += !run_test("shock_tube4.par") ;
fail_count += !run_test("shock_tube5.par") ;
fail_count += !run_test("shock_tube6.par") ;
fail_count += !run_test("shock_tube7.par") ;

count += !run_test("soundwave_32.par") ;
count += !run_test("soundwave_64.par") ;
count += !run_test("soundwave_128.par") ;
count += !run_test("soundwave_256.par") ;
count += !run_test("soundwave_512.par") ;
fail_count += !run_test("soundwave_32.par") ;
fail_count += !run_test("soundwave_64.par") ;
fail_count += !run_test("soundwave_128.par") ;
fail_count += !run_test("soundwave_256.par") ;
fail_count += !run_test("soundwave_512.par") ;

count += !run_test("friction_2spc.par", "friction_2spc.spc") ;
count += !run_test("friction_2spc_an.par", "friction_2spc.spc") ;
count += !run_test("friction_2spc_phys.par", "friction_2spc.spc") ;
count += !run_test("friction_2spc_phys_an.par", "friction_2spc.spc") ;
count += !run_test("friction_6spc.par", "friction_6spc.spc") ;
fail_count += !run_test("friction_2spc.par", "friction_2spc.spc") ;
fail_count += !run_test("friction_2spc_an.par", "friction_2spc.spc") ;
fail_count += !run_test("friction_2spc_phys.par", "friction_2spc.spc") ;
fail_count += !run_test("friction_2spc_phys_an.par", "friction_2spc.spc") ;
fail_count += !run_test("friction_6spc.par", "friction_6spc.spc") ;

count += !run_test("collheat_2spc.par", "collheat_2spc.spc") ;
count += !run_test("collheat_2spc_rad.par", "collheat_2spc.spc") ;
fail_count += !run_test("collheat_2spc.par", "collheat_2spc.spc") ;
fail_count += !run_test("collheat_2spc_rad.par", "collheat_2spc.spc") ;

count += !run_test("dustywave_nonstiff.par") ;
count += !run_test("dustywave_stiff.par") ;
fail_count += !run_test("dustywave_nonstiff.par") ;
fail_count += !run_test("dustywave_stiff.par") ;

count += !run_test("dusty_shock.par") ;
fail_count += !run_test("dusty_shock.par") ;

count += !run_test("planet_cartesian.par") ;
count += !run_test("planet_spherical.par") ;

count += !run_test("irradiation.par", "irradiation.spc") ;
fail_count += !run_test("planet_cartesian.par") ;
fail_count += !run_test("planet_spherical.par") ;

fail_count += !run_test("irradiation1.par", "irradiation.spc") ;
fail_count += !run_test("irradiation2.par", "irradiation.spc") ;

return count ;
return fail_count ;
}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand Down
Loading

0 comments on commit a65382b

Please sign in to comment.