[Richards] Enable/Disable flux check depending on FE order

parent d63a7681
......@@ -22,10 +22,12 @@ def evaluate(iniinfo,runtime):
initial_head = float(iniinfo["_ode.headLower"])
extensions = list(map(float,iniinfo["grid.extensions"].split()))
do_flux_rt = iniinfo["richards.fluxReconstruction.enable"]
fe_order = int(iniinfo["richards.numerics.FEorder"])
do_flux_rt = iniinfo["richards.fluxReconstruction.enable"]
do_flux_grad = fe_order > 0
if do_flux_rt:
flux_rt_key = "flux_RT" + str(max(0,fe_order-1))
# get parameter data
param_data = read_yml(iniinfo["richards.parameters.file"])
......@@ -117,19 +119,26 @@ def evaluate(iniinfo,runtime):
plt.savefig("{}/head_residual.png".format(iniinfo["richards.output.outputPath"]))
plt.close()
flux = data["flux"][iy,:]
res_flux = flux - np.array([0.,influx,0.])[np.newaxis,:]
l2_flux = [np.sqrt(average(np.square(r[iy_inv]))) for r in res_flux.T]
print("flux L2 errors: {:.2e}, {:.2e}, {:.2e}".format(*l2_flux))
# Check matric head
head_tol = 1E-5 if not "_ode.head_abstol" in iniinfo else float(iniinfo["_ode.head_abstol"])
passed = bool(l2_head < head_tol)
plt.figure()
plt.plot(y,res_flux)
plt.savefig("{}/flux_residual.png".format(iniinfo["richards.output.outputPath"]))
plt.close()
# Check flux from gradient
if do_flux_grad:
flux = data["flux"][iy,:]
res_flux = flux - np.array([0.,influx,0.])[np.newaxis,:]
l2_flux = [np.sqrt(average(np.square(r[iy_inv]))) for r in res_flux.T]
print("flux L2 errors: {:.2e}, {:.2e}, {:.2e}".format(*l2_flux))
plt.figure()
plt.plot(y,res_flux)
plt.savefig("{}/flux_residual.png".format(iniinfo["richards.output.outputPath"]))
plt.close()
head_tol = 1E-5 if not "_ode.head_abstol" in iniinfo else float(iniinfo["_ode.head_abstol"])
flux_tol = abs(influx) * 1e-5 if not "_ode.flux_abstol" in iniinfo else float(iniinfo["_ode.flux_abstol"])
flux_tol = abs(influx) * 1e-5 if not "_ode.flux_abstol" in iniinfo else float(iniinfo["_ode.flux_abstol"])
passed = passed and all([l < flux_tol for l in l2_flux])
# Check reconstructed flux
if do_flux_rt:
flux_rt = data[flux_rt_key][iy,:]
res_flux_rt = flux_rt - np.array([0.,influx,0.])[np.newaxis,:]
......@@ -143,9 +152,6 @@ def evaluate(iniinfo,runtime):
flux_rt_tol = 1e-13 if not "_ode.flux_rt_abstol" in iniinfo else float(iniinfo["_ode.flux_rt_abstol"])
return bool(l2_head < head_tol) \
and all([l < flux_tol for l in l2_flux]) \
and all([l < flux_rt_tol for l in l2_flux_rt])
passed = passed and all([l < flux_rt_tol for l in l2_flux_rt])
# no flux reconstruction
return bool(l2_head < head_tol) and all([l < flux_tol for l in l2_flux])
return passed
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment