diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 7b9b6698cea83c4b215120bf666db5b2655e04bc..455b2e72fbfb4756b598359be234f66c0efa5860 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -162,9 +162,9 @@ class LDDsimulation(object):
df.info(df.parameters, True)
- self.solver_type_is_Kryov = True
+ self.solver_type_is_Krylov = True
### Define the linear solver to be used.
- if self.solver_type_is_Kryov:
+ if self.solver_type_is_Krylov:
self.solver = 'bicgstab'#'superlu' #'gmres'#'bicgstab' # biconjugate gradient stabilized method
self.preconditioner = 'jacobi' #'hypre_amg' 'default' #'ilu'#'hypre_amg' # incomplete LU factorization
# dictionary of solver parametrs. This is passed to self._init_subdomains,
@@ -230,6 +230,8 @@ class LDDsimulation(object):
""" set parameters of an instance of class LDDsimulation"""
self.use_case = use_case
+ if include_gravity:
+ self.use_case = self.use_case + "-with-gravity"
self.output_dir = output_dir
self.subdomain_def_points = subdomain_def_points
self.isRichards = isRichards
@@ -352,6 +354,9 @@ class LDDsimulation(object):
self.t += self.timestep_size
# write solutions to files.
for subdom_ind, subdomain in self.subdomain.items():
+ # for phase in subdomain.has_phases:
+ # print(f"calculated pressure of {phase}phase:")
+ # print(subdomain.pressure[phase].vector()[:])
subdomain.write_solution_to_xdmf(#
file = self.solution_file[subdom_ind], #
time = self.t,#
@@ -361,11 +366,29 @@ class LDDsimulation(object):
# absolute error.
if subdomain.pressure_exact is not None:
relative_L2_errornorm = dict()
+ pressure_exact = dict()
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]
pa_exact.t = self.t
norm_pa_exact = df.norm(pa_exact, norm_type='L2', mesh=subdomain.mesh)
error_calculated = df.errornorm(pa_exact, subdomain.pressure[phase], 'L2', degree_rise=6)
+ pressure_exact.update(
+ {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
+ )
+ absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
+ pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
+ abs_diff_tmp = np.fabs(pressure_difference)
+ absolute_difference.vector()[:] = abs_diff_tmp
+ if phase == "wetting":
+ data_set_name="abs_diff_w"
+ else:
+ data_set_name="abs_diff_nw"
+ self.write_function_to_xdmf(#
+ function=absolute_difference,
+ file=self.solution_file[subdom_ind], #
+ time=self.t,#
+ data_set_name=data_set_name
+ )
if norm_pa_exact > self.tol:
relative_L2_errornorm.update(
{phase: error_calculated/norm_pa_exact}
@@ -385,7 +408,7 @@ class LDDsimulation(object):
self.timestep_num += 1
print("LDD simulation use case: {}".format(self.use_case))
- df.info("Finished")
+ df.info("Finished calculation")
df.info("Elapsed time:"+str(timer.elapsed()[0]))
timer.stop()
# r_cons.close()
@@ -397,7 +420,9 @@ class LDDsimulation(object):
# timefile.write(df.timings)
df.dump_timings_to_xml(self.output_dir+"timings.xml", df.TimingClear.keep)
-
+ # df.info(colored("start post processing calculations ...\n", "yellow"))
+ # self.post_processing()
+ # df.info(colored("finished post processing calculations \nAll right. I'm Done.", "green"))
def LDDsolver(self, time: float = None, #
debug: bool = False, #
@@ -695,6 +720,7 @@ class LDDsimulation(object):
# set up solution file names for each subdomain.
self.output_filename_parameter_part = dict()
self.solution_file = dict()
+ self.postprocessed_solution_file = dict()
for subdom_ind, subdomain in self.subdomain.items():
tau = subdomain.timestep_size
L = subdomain.L
@@ -707,7 +733,6 @@ class LDDsimulation(object):
+ "_gravitiy_" + "{}".format(self.include_gravity)
name3 = "_lambda" + "{}".format(Lam) + "_Lw" + "{}".format(L["wetting"])
filename_param_part = name1 + name2 + name3
- filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
else:
Lam_w = subdomain.lambda_param['wetting']
Lam_nw = subdomain.lambda_param['nonwetting']
@@ -719,7 +744,9 @@ class LDDsimulation(object):
+"_lambda_nw" + "{}".format(Lam_nw)\
+"_Lw" + "{}".format(L["wetting"])\
+"_Lnw" + "{}".format(L["nonwetting"])
- filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
+
+ filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
+ post_filename = self.output_dir + filename_param_part + "_solution_post" +".xdmf"
# create solution files and populate dictionary.
self.output_filename_parameter_part.update(
{subdom_ind: filename_param_part}
@@ -727,6 +754,10 @@ class LDDsimulation(object):
self.solution_file.update( #
{subdom_ind: SolutionFile(self._mpi_communicator, filename)}
)
+ # self.postprocessed_solution_file.update(
+ # {subdom_ind: SolutionFile(self._mpi_communicator, post_filename)}
+ # )
+
def write_exact_solution_to_xdmf(self):
"""
@@ -736,14 +767,26 @@ class LDDsimulation(object):
for subdom_ind, subdomain in self.subdomain.items():
t = copy.copy(self.starttime)
file = self.solution_file[subdom_ind]
- for timestep in range(0,self.number_of_timesteps+1):
+
+ for timestep in range(self.number_of_timesteps+1):
+ exact_pressure = dict()
for phase in subdomain.has_phases:
pa_exact = subdomain.pressure_exact[phase]
pa_exact.t = t
- tmp_exact_pressure = df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])
- tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
- file.write(tmp_exact_pressure, t)
-
+ exact_pressure.update(
+ {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])}
+ )
+ exact_pressure[phase].rename("exact_pressure_"+"{}".format(phase), "exact_pressure_"+"{}".format(phase))
+ file.write(exact_pressure[phase], t)
+
+ exact_capillary_pressure = df.Function(subdomain.function_space["pressure"]['wetting'])
+ if self.isRichards:
+ exact_capillary_pressure.assign(-exact_pressure['wetting'])
+ else:
+ pc_temp = exact_pressure["nonwetting"].vector()[:] - exact_pressure["wetting"].vector()[:]
+ exact_capillary_pressure.vector()[:] = pc_temp
+ exact_capillary_pressure.rename("pc_exact", "pc_exact")
+ file.write(exact_capillary_pressure, t)
t += self.timestep_size
@@ -842,6 +885,28 @@ class LDDsimulation(object):
else:
self.errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)
+
+ def write_function_to_xdmf(self,
+ function: df.Function,
+ file: tp.Type[SolutionFile], #
+ data_set_name: str,
+ time: float = None,#
+ ):
+ """
+ Write function to disk file file. It is expected, that file is a df.XDMFFile
+
+ PARAMETERS
+ function df.Function dolfin function containing the
+ information to write to the file.
+ file str File name of df.XDMFFile to write to
+ data_set_name str name of the data set
+ time float time at which to write the solution.
+
+ """
+ function.rename("{}".format(data_set_name), "{}".format(data_set_name))
+ file.write(function, time)
+
+
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
mesh_resolution: float = None, #
@@ -1044,7 +1109,7 @@ class LDDsimulation(object):
)
# setup the linear solvers and set the solver parameters.
# df.LinearSolver(self.solver)
- if self.solver_type_is_Kryov:
+ if self.solver_type_is_Krylov:
self.subdomain[subdom_num].linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
else:
self.subdomain[subdom_num].linear_solver = df.LUSolver()
@@ -1054,13 +1119,13 @@ class LDDsimulation(object):
# assign solver parameters set at the beginning of the LDDsimulation class.
solver_param = self.subdomain[subdom_num].linear_solver.parameters
- df.info(solver_param, True)
+ df.info(solver_param, False)
for parameter_name in self.solver_parameters.keys():
solver_param[parameter_name] = self.solver_parameters[parameter_name]
- # ## print out set solver parameters
- for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
- print(f"parameter: {parameter} = {value}")
+ # # ## print out set solver parameters
+ # for parameter, value in self.subdomain[subdom_num].linear_solver.parameters.items():
+ # print(f"parameter: {parameter} = {value}")
if self.write2file['meshes_and_markers']:
filepath = self.output_dir+f"interface_marker_subdomain{subdom_num}.pvd"
@@ -1129,6 +1194,22 @@ class LDDsimulation(object):
subdomain.pressure[phase].assign(pa0_interpolated)
subdomain.pressure_prev_iter[phase].assign(pa0_interpolated)
subdomain.pressure_prev_timestep[phase].assign(pa0_interpolated)
+ # write zeros to the abs_diff variables
+
+ absolute_difference = df.Function(subdomain.function_space["pressure"][phase])
+ # pressure_difference = pressure_exact[phase].vector()[:] - subdomain.pressure[phase].vector()[:]
+ # abs_diff_tmp = np.fabs(pressure_difference)
+ # absolute_difference.vector()[:] = abs_diff_tmp
+ if phase == "wetting":
+ data_set_name="abs_diff_w"
+ else:
+ data_set_name="abs_diff_nw"
+ self.write_function_to_xdmf(#
+ function=absolute_difference,
+ file=self.solution_file[num], #
+ time=self.starttime,#
+ data_set_name=data_set_name
+ )
# populate interface dictionaries with pressure values. The calc_gli_term
# of each subdomain reads from the interface.pressure_values_prev
@@ -1229,6 +1310,65 @@ class LDDsimulation(object):
)
+ # def post_processing(self):
+ # """post processing of the simulation.
+ # calculate
+ # - pc_exact and pc_num
+ # - absolute differences to exact solution if present.
+ # - relative errors to exact solution if present
+ # """
+ # for subdom_ind, subdomain in self.subdomain.items():
+ # self._mesh = df.Mesh()
+ # solution_file = self.solution_file[subdom_ind]
+ # post_solution_file = self.postprocessed_solution_file[subdom_ind]
+ # pressure = dict()
+ # # the internal time series counter in df.XDMFFile starts at 0 and
+ # # refers to the first saved value. This corresponds to the initial
+ # # values. We then calculated self.number_of_timesteps many timesteps.
+ # t = self.starttime
+ # for internal_timestep in range(self.number_of_timesteps+1):
+ # for phase in subdomain.has_phases:
+ # pressure.update(
+ # {phase: df.Function(subdomain.function_space["pressure"][phase])}
+ # )
+ # phase_pressure_string = "pressure_"+"{}".format(phase)
+ # solution_file.read_checkpoint(
+ # u=pressure[phase],
+ # name=phase_pressure_string,
+ # counter=internal_timestep
+ # )
+ # pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
+ # post_solution_file.write(pressure[phase], t)
+ # t += self.timestep_size
+ # print(f"read pressure of {phase}phase:")
+ # print(pressure[phase].vector()[:])
+ # solution_file.close()
+ # # # if we have an exact solution, write out the |u -uh|_L2 to see the
+ # # # absolute error.
+ # # if subdomain.pressure_exact is not None:
+ # # relative_L2_errornorm = dict()
+ # # for phase in subdomain.has_phases:
+ # # pa_exact = subdomain.pressure_exact[phase]
+ # # pa_exact.t = self.t
+ # # norm_pa_exact = df.norm(pa_exact, norm_type='L2', mesh=subdomain.mesh)
+ # # error_calculated = df.errornorm(pa_exact, subdomain.pressure[phase], 'L2', degree_rise=6)
+ # # if norm_pa_exact > self.tol:
+ # # relative_L2_errornorm.update(
+ # # {phase: error_calculated/norm_pa_exact}
+ # # )
+ # # else:
+ # # relative_L2_errornorm.update(
+ # # {phase: error_calculated}
+ # # )
+ # # errornorm_filename = self.output_dir+self.output_filename_parameter_part[subdom_ind]+\
+ # # "_L2_errornorms_over_time" +".csv"
+ # # self.write_errornorms_to_csv(
+ # # filename = errornorm_filename, #
+ # # subdomain_index = subdom_ind,
+ # # errors = relative_L2_errornorm,
+ # # )
+
+
def _init_exact_solution_expression(self):
""" initialise dolfin expressions for exact solution and populate
self.exact_solution_expr dictionary. This is used to calculate errornorms
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index dedbec97db1323c3f16af2a8d03f322eb599039f..9204d77c5ea7416fd70b5004813ebc739676a883 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -363,9 +363,21 @@ class Interface(BoundaryPart):
for mesh_entity in mesh_entity_iterator:
# this is the facet index relative to the Mesh
# that interface_marker is defined on.
- facet_index = mesh_entity.index()
+ interface_facet_index = mesh_entity.index()
+ for interface_facet in df.facets(mesh_entity):
+ # we only need the x,y coordinates. vertex.point()[:] returns
+ # the array and cannot be sliced. the result can be indexed, hence,
+ # vertex.point()[:][0:2]
+ # facet_midpoint = interface_facet.midpoint()[:][0:2]
+ # print(f"facet: {interface_facet}")
+ # print("facet coordinates: ")
+ facet_points = []
+ for vertex in df.vertices(interface_facet):
+ # print(f"vertex: {vertex} with coordinates {vertex.point()[:][0:2]}")
+ facet_points.append(vertex.point())
+
- dofs_and_coords.update({facet_index: dict()})
+ dofs_and_coords.update({interface_facet_index: dict()})
# because mesh_entity doesn't have a method to access coordinates we
# need to cast it into an acutal facet object.
for cell in df.cells(mesh_entity):
@@ -374,6 +386,8 @@ class Interface(BoundaryPart):
# returns an iterator object.
facet_cell = cell
facet_cell_index = facet_cell.index()
+ # cell_facets_entity = df.facets(cell)
+ # cell_coordinates = cell.get_vertex_coordinates()
# the cell is a triangle and for each dof numbered locally for that
# cell, we have the coordinates in a numpy array. The numbering is the
@@ -389,9 +403,9 @@ class Interface(BoundaryPart):
# cell is a triangle and only has one facet that belongs to
# the interface. we want to store only dofs that lie on that particular
# facet
- if self._is_on_boundary_part(dof_coord):
+ if self._is_on_line_segment(facet_points[0], facet_points[1], dof_coord):
global_dof_ind = facet_cell_dofmap[local_dof_ind]
- dofs_and_coords[facet_index].update({global_dof_ind: dof_coord})
+ dofs_and_coords[interface_facet_index].update({global_dof_ind: dof_coord})
if debug:
@@ -484,10 +498,11 @@ class Interface(BoundaryPart):
for local_facet_index, local_facet_vertex in subdom_facet2coord.items():
for global_facet_index, global_facet_vertex in global_facet2coord.items():
# vertices are not necessarily in the same order
- vertices_are_equal = np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
- vertices_are_equal += np.array_equal(local_facet_vertex[0], global_facet_vertex[1])
- vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[0])
- vertices_are_equal += np.array_equal(local_facet_vertex[1], global_facet_vertex[1])
+ # np.array_equal(local_facet_vertex[0], global_facet_vertex[0])
+ vertices_are_equal = np.allclose(local_facet_vertex[0], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
+ vertices_are_equal += np.allclose(local_facet_vertex[0], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
+ vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[0], rtol=1e-10, atol=1e-14)
+ vertices_are_equal += np.allclose(local_facet_vertex[1], global_facet_vertex[1], rtol=1e-10, atol=1e-14)
both_vertices_are_equal = (vertices_are_equal == 2)
if both_vertices_are_equal:
self.local_to_global_facet_map[subdomain_index].update(
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 380c3414ebd158473e391d45a738c1397e8bbc90..0b11785b48f362ca4e85c24eb0f41666b7aec23e 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -388,7 +388,7 @@ class DomainPatch(df.SubDomain):
# the non wetting phase has a + before instead of a - before the bracket
rhs2 = (porosity*(S(pc_prev_iter) - S(pc_prev_timestep))*phi_a)*dx
if self.include_gravity:
- # z_a included the minus sign already, see self._calc_gravity_expressions
+ # z_a included the minus sign, see self._calc_gravity_expressions
rhs_gravity = dt/mu_a*df.dot(ka(1 - S(pc_prev_iter))*df.grad(z_a), df.grad(phi_a))*dx
else:
raise RuntimeError('missmatch for input parameter phase.',
@@ -398,14 +398,14 @@ class DomainPatch(df.SubDomain):
form = form1 + form2 + sum(interface_forms)
if self.include_gravity:
# z_a included the minus sign already, see self._calc_gravity_expressions
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) - rhs_gravity
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
else:
rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
else:
form = form1 + form2
if self.include_gravity:
# z_a included the minus sign already, see self._calc_gravity_expressions
- rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - rhs_gravity
+ rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx + rhs_gravity
else:
rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx
form_and_rhs = {#
@@ -468,15 +468,15 @@ class DomainPatch(df.SubDomain):
# to calculate a gl0 term.
# TODO: add intrinsic permeabilty!!!!
if self.include_gravity:
- # z_a included the minus sign already, see self._calc_gravity_expressions
- flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev + z_a)
+ # z_a included the minus sign, see self._calc_gravity_expressions
+ flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev - z_a)
else:
flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
else:
# phase == 'nonwetting'
# we need anothre flux in this case
if self.include_gravity:
- flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev + z_a)
+ flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
else:
flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
@@ -603,20 +603,42 @@ class DomainPatch(df.SubDomain):
gli_prev = interface.gli_term_prev[neighbour][phase][global_facet_index]
p_prev = pressure_prev[global_facet_index]
gli_save_dict = gli_save_dictionary[global_facet_index]
+ # print(f"\non subdomain{self.subdomain_index}")
+ # print(f"on facet with glob_index: {global_facet_index} we have gli_dofs")
+ #
for gli_dof_index, gli_dof_coord in gli_facet_dofs2coordinates.items():
# find the right previous pressure dof of the interface
# dictionary to use for the calculation.
+ # print(f"coordinates of gli_dof with index {gli_dof_index} I try to match is: {gli_dof_coord}" )
p_previous_dof = None
for p_dof_coord_tupel, p_prev_dof in p_prev.items():
p_coord = np.array([p_dof_coord_tupel[0], p_dof_coord_tupel[1]])
- if np.array_equal(gli_dof_coord, p_coord):
+ # print(f"coordinates of p_dof I test this against: {p_coord}")
+ a_close_b = np.allclose(gli_dof_coord, p_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(p_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ # print(f"type(gli_dof_coord) = {type(gli_dof_coord)}")
+ # print(f"np.allclose({gli_dof_coord}, {p_coord}, rtol=1e-10, atol=1e-14): {a_close_b }")
+ # print(f"np.allclose({p_coord}, {gli_dof_coord}, rtol=1e-10, atol=1e-14): {b_close_a}")
+ if a_close_b and b_close_a:
p_previous_dof = p_prev_dof
+ print("\n")
+ if p_previous_dof is None:
+ raise RuntimeError(f"didn't find p_previous_dof corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
# same with the gli_previous_dof
gli_previous_dof = None
for gli_prev_coord_tupel, gli_prev_dof in gli_prev.items():
gli_prev_coord = np.array([gli_prev_coord_tupel[0], gli_prev_coord_tupel[1]])
- if np.array_equal(gli_dof_coord, gli_prev_coord):
+ # due to the warning in the documentation of np.allclose
+ # that np.allclose is not symetric, we test if both
+ # np.allclose(a,b) and np.allclose(b,a) are true.
+ a_close_b = np.allclose(gli_dof_coord, gli_prev_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(gli_prev_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
gli_previous_dof = gli_prev_dof
+ if gli_previous_dof is None:
+ raise RuntimeError(f"didn't find gli_previous_dof corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
gli_value = -2*Lambda*p_previous_dof - gli_previous_dof
gli.vector()[gli_dof_index] = gli_value
@@ -871,25 +893,42 @@ class DomainPatch(df.SubDomain):
# dof wise.
for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind].update(
- {gli_dof_index: dict()}
+ {gli_dof_index: {"flux_x": None,
+ "flux_x": None,
+ "pressure": None}}
)
for flux_x_dof_ind, dof_coord in flux_x_dofs_along_facet.items():
- if np.array_equal(dof_coord, gli_dof_coord):
+ a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
{"flux_x": flux_x_dof_ind}
)
+ if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_x"] is None:
+ raise RuntimeError(f"didn't find flux_x_dof_ind corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
for flux_y_dof_ind, dof_coord in flux_y_dofs_along_facet.items():
- if np.array_equal(dof_coord, gli_dof_coord):
+ a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
{"flux_y": flux_y_dof_ind}
)
+ if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["flux_y"] is None:
+ raise RuntimeError(f"didn't find flux_y_dof_ind corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
for p_dof_ind, dof_coord in p_dofs_along_facet.items():
- if np.array_equal(dof_coord, gli_dof_coord):
+ a_close_b = np.allclose(gli_dof_coord, dof_coord, rtol=1e-10, atol=1e-14)
+ b_close_a = np.allclose(dof_coord, gli_dof_coord, rtol=1e-10, atol=1e-14)
+ if a_close_b and b_close_a:
self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index].update(
{"pressure": p_dof_ind}
)
+ if self.interface_corresponding_dof_index[interf_ind][phase][local_facet_ind][gli_dof_index]["pressure"] is None:
+ raise RuntimeError(f"didn't find p_dof_ind corresponding to dict key ({gli_dof_coord})",
+ "something is wrong")
def _init_interface_values(self):
@@ -938,17 +977,24 @@ class DomainPatch(df.SubDomain):
write_iter_for_fixed_time bool flag to toggle wether pairs of
(solution, iteration) for fixed t should be
written instead of (solution, t)
- subsequent_errors df.Funtion FEM Function passed by the solver
- containing self.pressure - self.pressure_prev_iter
"""
t = time
- for phase in self.has_phases:
- if not write_iter_for_fixed_time:
+ if not write_iter_for_fixed_time:
+ for phase in self.has_phases:
self.pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
file.write(self.pressure[phase], t)
+ capillary_pressure = df.Function(self.function_space["pressure"]["wetting"])
+ if self.isRichards:
+ capillary_pressure.assign(-self.pressure["wetting"])
else:
+ pc_temp = self.pressure["nonwetting"].vector()[:] - self.pressure["wetting"].vector()[:]
+ capillary_pressure.vector()[:] = pc_temp
+ capillary_pressure.rename("pc_num", "pc_num")
+ file.write(capillary_pressure, t)
+ else:
+ for phase in self.has_phases:
i = self.iteration_number
self.pressure[phase].rename("pressure_"+"{}".format(phase), #
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
@@ -961,9 +1007,5 @@ class DomainPatch(df.SubDomain):
"pressure(t="+"{}".format(t)+")_"+"{}".format(phase)
)
file.write(error, i)
- # self.gli_function[phase].vector().set_local(self.gli_assembled[phase])
- # self.gli_function[phase].rename("gli_function_"+"{}".format(phase),#
- # "all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase))
- # file.write(self.gli_function[phase], i)
# END OF CLASS DomainPatch
diff --git a/LDDsimulation/helpers.py b/LDDsimulation/helpers.py
index 803031f9ec1f51ff17d62d1b37172ed63ff372db..765630c9c623a3b85ec38414fdd65e44a7eab9f4 100644
--- a/LDDsimulation/helpers.py
+++ b/LDDsimulation/helpers.py
@@ -111,7 +111,7 @@ def generate_exact_solution_expressions(
dxdxflux = -1/mu*dka(1-S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(1-S)
# y part of div(flux) for nonwetting
if include_gravity:
- dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(1-S)
+ dydyflux = -1/mu*dka(1-S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(1-S)
else:
dydyflux = -1/mu*dka(1-S)*dS*dypc*dypa + 1/mu*dydypa*ka(1-S)
else:
@@ -119,7 +119,7 @@ def generate_exact_solution_expressions(
dxdxflux = 1/mu*dka(S)*dS*dxpc*dxpa + 1/mu*dxdxpa*ka(S)
# y part of div(flux) for wetting
if include_gravity:
- dydyflux = 1/mu*dka(S)*dS*dypc*(dypa - rho*g) + 1/mu*dydypa*ka(S)
+ dydyflux = 1/mu*dka(S)*dS*dypc*(dypa + rho*g) + 1/mu*dydypa*ka(S)
else:
dydyflux = 1/mu*dka(S)*dS*dypc*(dypa) + 1/mu*dydypa*ka(S)