Commit 0d810dba authored by Dion Haefner's avatar Dion Haefner

added evaporation ode, parallel correlation and interpolator ode

parent bc4d4db4
......@@ -35,7 +35,7 @@ def evaluate(iniinfo,runtime):
return False
for key in f1_data.keys():
if not np.array_equal(f1_data[key],f2_data[key]):
if not np.allclose(f1_data[key],f2_data[key]):
print("Data does not match in key {}".format(key))
return False
......
......@@ -50,18 +50,30 @@ def evaluate(iniinfo,runtime):
t0 = 0 # y-coordinate of grid bottom
tmax = extensions[1] # vertical extensions of the grid
if parfield_data["raw_field"].shape == (2,1):
edge = .5*tmax
if parfield_data["raw_field"].size == 1:
par = {k: parfield_data[k][0] for k in parfield_data.keys()}
parfun = lambda y: par
elif parfield_data["raw_field"].shape == (2,1):
par = [{k: parfield_data[k][l,0] for k in parfield_data.keys()} for l in range(2)]
elif parfield_data["raw_field"].size == 1:
edge = np.inf # only 1 medium -> push edge to infinity
par = [{k: parfield_data[k][0] for k in parfield_data.keys()},None]
parfun = lambda y: par[0] if y <= 0.5*extensions[1] else par[1]
elif parfield_data["raw_field"].shape == (2,2) and iniinfo["parameters.interpolation"] == "linear":
assert all(np.allclose(v[:,0],v[:,1]) for v in parfield_data.values())
par = [{k: parfield_data[k][l,0] for k in parfield_data.keys()} for l in range(2)]
def parfun(y):
y_norm = y / extensions[1]
if y_norm < 0.25:
return par[0]
elif y_norm < 0.75:
y_scale = 2 * (y_norm - 0.25)
return {k: par[0][k] * (1-y_scale) + par[1][k] * y_scale for k in parfield_data.keys()}
else:
return par[1]
else:
raise ValueError("Parameter field must be homogeneous or layered vertically (1 or 2 cells)")
raise ValueError("Parameter field must be homogeneous or layered vertically")
print("Operating on file: {}".format(vtk_file))
grid = vtk_file.read(["head","flux"])
grid = vtk_file.read(["head","flux","K_0"])
cells = grid.cellCenters()
y = cells[...,1].flatten()
tri, connectivity = grid.triangulation()
......@@ -70,10 +82,18 @@ def evaluate(iniinfo,runtime):
y = y[iy]
head = data["head"].flatten()[iy]
k0_dorie = data["K_0"].flatten()[iy]
k0_calc = [parfun(yi)["k0"] for yi in y]
if not np.allclose(k0_dorie,k0_calc):
print("Encountered unexpected parameter values")
return False
del k0_dorie
del k0_calc
r = integrate.ode(richards_equation.richards)
r.set_integrator('dopri5')
r.set_integrator("dopri5", atol=1e-10, rtol=1e-10)
r.set_initial_value(initial_head, 0)
r.set_f_params(influx,par[0],par[1],edge)
r.set_f_params(influx,parfun)
r_y = []
for t in y:
......
......@@ -17,9 +17,9 @@ def k(h,k0,tau,alpha,n,**kwargs):
"""
return k0*np.power(1+np.power(alpha*h,n),-tau*(1-1/n))*np.power(1-np.power(1-1/(1+np.power(alpha*h,n)),1-1/n),2)
def richards(y,h,jw,param1,param2,edge):
def richards(y,h,jw,param):
"""
Returns :math:`du/dy` as given by the Richards equation in a 1-dimensional, layered medium with constant inflow:
Returns :math:`du/dy` as given by the Richards equation in a 1-dimensional medium with constant inflow:
.. math::
......@@ -33,11 +33,7 @@ def richards(y,h,jw,param1,param2,edge):
:param float edge: Position of the edge dividing the two sets of parameters
"""
if (y<=edge):
a=k(h,**param1)
else:
a=k(h,**param2)
return -1-jw/a
return -1-jw/k(h,**param(y))
def van_genuchten(h,alpha,n,theta_r,theta_s):
"""
......
......@@ -12,23 +12,25 @@ function(dorie_add_system_test_dependency test1 test2)
endfunction()
# dorie run
#dorie_add_system_test(dorie ode.mini)
#dorie_add_system_test(dorie ode_higherorder.mini)
#dorie_add_system_test(dorie muphi.mini)
#dorie_add_system_test_dependency(dorie_muphi dorie-pfg_muphi_pfg)
#dorie_add_system_test(dorie reference_2d.mini)
#dorie_add_system_test(dorie reference_3d.mini)
#dorie_add_system_test(dorie reference_evaporation.mini)
dorie_add_system_test(dorie ode.mini)
dorie_add_system_test(dorie ode_higherorder.mini)
dorie_add_system_test(dorie ode_evaporation.mini)
dorie_add_system_test(dorie ode_linear_interpolator.mini)
dorie_add_system_test(dorie muphi.mini)
dorie_add_system_test_dependency(dorie_muphi dorie-pfg_muphi_pfg)
dorie_add_system_test(dorie reference_2d.mini)
dorie_add_system_test(dorie reference_3d.mini)
dorie_add_system_test(dorie reference_evaporation.mini)
# dorie pfg
#dorie_add_system_test(dorie-pfg parfield.mini)
dorie_add_system_test(dorie-pfg parfield_parallel.mini)
#dorie_add_system_test(dorie-pfg parfield_muphi.mini)
#dorie_add_system_test(dorie-pfg correlation.mini)
dorie_add_system_test(dorie-pfg parfield.mini)
dorie_add_system_test(dorie-pfg parfield_muphi.mini)
dorie_add_system_test(dorie-pfg correlation.mini)
dorie_add_system_test(dorie-pfg correlation_parallel.mini)
# dorie plot
#dorie_add_system_test(dorie plot.mini)
#dorie_add_system_test_dependency(dorie_plot dorie_reference_2d_0000)
dorie_add_system_test(dorie plot.mini)
dorie_add_system_test_dependency(dorie_plot dorie_reference_2d_0000)
# dorie create
#dorie_add_system_test(dorie create.mini)
dorie_add_system_test(dorie create.mini)
......@@ -2,7 +2,5 @@ spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 3
0 neumann -5.55e-6 dirichlet 0
1e5 neumann 0 dirichlet 0
1e5 evaporation 1e-7 -10 dirichlet 0
number_BC_change_times 1
0 evaporation 2e-7 -100 dirichlet 0
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 3
0 neumann -5.55e-6 dirichlet 0
1e5 neumann 0 dirichlet 0
1e5 evaporation 1e-7 -10 dirichlet 0
include ${CMAKE_BINARY_DIR}/doc/default_files/parfield.ini
__name = pfg_parallel
__name = correlation_parallel
_test_command = pfg
_test_command_options = --parallel -m=--allow-run-as-root
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
......@@ -19,7 +19,7 @@ dimensions = 2, 3 | expand dim
extensions = 1 1, 1 1 1 | expand dim
[generator.fft]
outputPath = fft_fields_parallel | unique
outputPath = correlation_parallel | unique
N = 2048 2048, 128 128 128 | expand dim
variance = 0.2
correlationLengths = .004 .005, .05 .06 .08 | expand dim
......
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
__name = ode_evaporation
_test_command = run
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = ode
output.fileName = ode_evaporation | unique
output.outputPath = ode_evaporation | unique
output.verbose = 0
time.end = 1E7
time.maxTimestep = 1E7
time.startTimestep = 1E4
adaptivity.threshold = 1E-10
adaptivity.useAdaptivity = false, true | expand distribution
grid.gridType = rectangular, gmsh | expand gridtype
grid.cells = 1 100
grid.initialLevel = 3, 0 | expand distribution
grid.gridFile = "{_asset_path}/meshes/square.msh"
boundary.file = "{_asset_path}/bcs/evaporation_2d.dat"
parameters.arrayFile = "{_asset_path}/parfields/sand.h5", "{_asset_path}/parfields/silt.h5" | expand distribution
[_ode]
flux = 2e-7
head_abstol = 1E-5
flux_abstol = 1E-9
[NewtonParameters]
AbsoluteLimit = 1E-10
Reduction = 1E-10
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
__name = ode_linear_interpolator
_test_command = run
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_evaluation = ode
output.fileName = ode_linear_interpolator | unique
output.outputPath = ode_linear_interpolator | unique
output.verbose = 0
time.end = 1E7
time.maxTimestep = 1E7
time.startTimestep = 1E4
adaptivity.threshold = 1E-8
adaptivity.useAdaptivity = false, true | expand gridtype
grid.gridType = rectangular, gmsh | expand gridtype
grid.cells = 1 100
grid.initialLevel = 2, 0 | expand gridtype
grid.gridFile = "{_asset_path}/meshes/square.msh"
boundary.file = "{_asset_path}/bcs/infiltration_2d.dat"
parameters.arrayFile = "{_asset_path}/parfields/corners.h5"
parameters.interpolation = linear
[_ode]
flux = -3e-6
head_abstol = 1E-5
flux_abstol = 1E-8
[NewtonParameters]
AbsoluteLimit = 1E-8
Reduction = 1E-8
......@@ -23,6 +23,6 @@ grid.gridType = gmsh, rectangular | expand
grid.cells = 20 20
grid.extensions = 1 1
grid.gridFile = "{_asset_path}/meshes/square.msh"
boundary.file = "{_asset_path}/bcs/evaporation_2d.dat"
boundary.file = "{_asset_path}/bcs/infiltration_evaporation_2d.dat"
parameters.arrayFile = "{_asset_path}/parfields/fft_2d.h5"
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