Skip to content

Commit

Permalink
Merge branch 'fix/pass-tests' into feature/parthenon-bump
Browse files Browse the repository at this point in the history
  • Loading branch information
Ben Prather committed Oct 20, 2023
2 parents 1370143 + 9750787 commit cc725b0
Show file tree
Hide file tree
Showing 8 changed files with 481 additions and 18 deletions.
3 changes: 2 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,11 @@ logs/
*.rhdf
*.xdmf
*.hst
# Archival files
# Archival/test files
kharma_parsed_*.par
log_*.txt
bondi_analytic_*.txt
atmosphere_soln_*.txt

# Editor documents
.project
Expand Down
2 changes: 1 addition & 1 deletion pars/emhd/bondi_viscous.par
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,7 @@ check_inflow_outer_x1 = false
check_inflow_inner_x1 = true
# Force outflow bounds for EMHD vars
outflow_EMHD_inner_x1 = true
outflow_EMHD_outer_x1 = false
outflow_EMHD_outer_x1 = true

<debug>
verbose = 1
Expand Down
26 changes: 13 additions & 13 deletions tests/bondi_viscous/check.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,40 +43,40 @@
state = bondi.get_bondi_fluid_state(mdot, rc, gam, dump.grid)
state.params['eta'] = eta
state.params['tau'] = tau
dP_check = bondi.compute_dP(mdot, rc, gam, dump.grid, eta, tau)
state.cache['dP'] = dP_check

# compute dP
# compute dP either by adjusting dump to include higher-order terms,
# or the computed state to exclude them
if dump['emhd/higher_order_terms']:
print("Res: "+str(res)+"; higher order terms enabled")
Theta = (dump['gam'] - 1.) * dump['u'] / dump['rho']
nu_emhd = eta / dump['rho']
dP = dump['dP'] * np.sqrt(nu_emhd * dump['rho'] * Theta / tau)
# we're directly modifying the cache here. Inadvisable
dump.cache['dP'] = dump['dP'] * np.sqrt(eta * Theta / tau)
state.cache['dP'] = bondi.compute_dP(mdot, rc, gam, dump.grid, eta, tau, start=np.mean(dump['dP'][-1]))
else:
dP = dump['dP']
Theta = (dump['gam'] - 1.) * dump['u'] / dump['rho']
nu_emhd = eta / dump['rho']
dP_check /= np.sqrt(nu_emhd * dump['rho'] * Theta / tau)
state.cache['dP'] = bondi.compute_dP(mdot, rc, gam, dump.grid, eta, tau, start=np.mean(dump['dP'][-1])) / \
np.sqrt(eta * Theta / tau)

# Plot
for var in ['rho', 'u', 'B1']:
for var in ['rho', 'u', 'B1', 'dP']:
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(1,1,1)
pplt.plot_diff_xz(ax, dump, state, var)
fig.savefig("compare_{}_{}.png".format(var, res))
plt.close(fig)

radius = np.mean(dump.grid['r'], axis=(1,2))
plt.plot(radius, dP_check, label='dP ODE')
plt.plot(radius, np.mean(dP, axis=(1,2)), label='dP code')
plt.plot(radius, state['dP'], label='dP ODE')
plt.plot(radius, np.mean(dump['dP'], axis=(1,2)), label='dP code')
plt.plot(radius, np.mean(dump['dP'], axis=(1,2)) - state['dP'], label='dP diff')
plt.legend()
plt.savefig('compare_dP_{}.png'.format(res))
plt.savefig('compare_dP1d_{}.png'.format(res))
plt.close()

# compute L1 norm
L1[r,0] = np.mean(np.fabs(dump['rho'] - state['rho']))
L1[r,1] = np.mean(np.fabs(dump['u'] - state['u']))
L1[r,2] = np.mean(np.fabs(dP - dP_check))
L1[r,2] = np.mean(np.fabs(np.mean(dump['dP'], axis=(1,2)) - state['dP'])[2:])
L1[r,3] = np.mean(np.fabs(dump['B1'] - state['B1']))

# MEASURE CONVERGENCE
Expand Down
Loading

0 comments on commit cc725b0

Please sign in to comment.