[Flux Reconstruction] Add comparision of flux reconstruction values in ode test

parent e580e288
......@@ -22,6 +22,11 @@ def evaluate(iniinfo,runtime):
initial_head = float(iniinfo["richards.initial.headLower"])
extensions = list(map(float,iniinfo["grid.extensions"].split()))
do_flux_rt = iniinfo["richards.numerics.fluxReconstruction"]
fe_order = int(iniinfo["richards.numerics.FEorder"])
if do_flux_rt:
flux_rt_key = "flux_RT" + str(fe_order-1)
# get parameter data
param_data = read_yml(iniinfo["richards.parameters.file"])
mapping_file = iniinfo["grid.mapping.file"]
......@@ -62,7 +67,10 @@ def evaluate(iniinfo,runtime):
print("Operating on file: {}".format(vtk_file))
grid = vtk_file.read(["head","flux","K_0"])
variables = ["head","flux","K_0","flux"]
if do_flux_rt:
variables.append(flux_rt_key)
grid = vtk_file.read(variables)
cells = grid.cellCenters()
y = cells[...,1].flatten()
tri, connectivity = grid.triangulation()
......@@ -114,6 +122,12 @@ def evaluate(iniinfo,runtime):
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))
if do_flux_rt:
flux_rt = data[flux_rt_key][iy,:]
res_flux_rt = flux_rt - np.array([0.,influx,0.])[np.newaxis,:]
l2_flux_rt = [np.sqrt(average(np.square(r[iy_inv]))) for r in res_flux_rt.T]
print(flux_rt_key + " L2 errors: {:.2e}, {:.2e}, {:.2e}".format(*l2_flux_rt))
plt.figure()
plt.plot(y,res_flux)
plt.savefig("{}/flux_residual.png".format(iniinfo["richards.output.outputPath"]))
......
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