Skip to content
Snippets Groups Projects
Commit c86096d1 authored by David Seus's avatar David Seus
Browse files

I might have a working code

parent e86796ce
No related branches found
No related tags found
No related merge requests found
......@@ -81,7 +81,7 @@ class LDDsimulation(object):
## Private variables
# maximal number of L-iterations that the LDD solver uses.
self._max_iter_num = 5
self._max_iter_num = 100
# TODO rewrite this with regard to the mesh sizes
self.calc_tol = self.tol
# The number of subdomains are counted by self.init_meshes_and_markers()
......@@ -180,7 +180,7 @@ class LDDsimulation(object):
self._init_subdomains()
self._init_initial_values()
def LDDsolver(self, time: float):
def LDDsolver(self, time: float, debug: bool = False):
""" calulate the solution on all patches for the next timestep
The heart of the simulation. The function LDDsolver implements the LDD
......@@ -215,7 +215,7 @@ class LDDsimulation(object):
subdomain.null_all_interface_iteration_numbers()
# update the time of the source terms before starting the iteration.
# The sources and sinks are the same throughout the iteration.
self._eval_sources(interpolation_degree = 2,
self._eval_sources(interpolation_degree = 3,
time = time,
subdomain_index = ind)
self.update_DirichletBC_dictionary(time = time, subdomain_index = ind)
......@@ -228,8 +228,8 @@ class LDDsimulation(object):
subdomain.iteration_number = 0
for phase in subdomain.has_phases:
# set the solution of the previous timestep as new prev_timestep pressure
subdomain.pressure_prev_timestep[phase].vector().set_local(
subdomain.pressure[phase].vector().get_local()
subdomain.pressure_prev_timestep[phase].assign(
subdomain.pressure[phase]
)
### actual iteration starts here
......@@ -246,25 +246,32 @@ class LDDsimulation(object):
if not subdomain_calc_isdone[sd_index]:
# set subdomain iteration to current iteration
subdomain.iteration_number = iteration
if debug:
print(f"Entering iteration {iteration}")
# solve the problem on subdomain
self.Lsolver_step(subdomain_index = sd_index,#
debug = True,
)
# subsequent_iterations_err = self.calc_iteration_error(
# subdomain_index = sd_index,#
# error_type = "L2"
# )
# # check if error criterion has been met
# if subsequent_iterations_err < self.calc_tol:
# subdomain_calc_isdone[sd_index] = True
debug = debug,
)
subsequent_iter_error = dict()
for phase in subdomain.has_phases:
pa = subdomain.pressure[phase]
pa_prev_i = subdomain.pressure_prev_iter[phase]
subsequent_iter_error.update(
{phase: df.errornorm(pa, pa_prev_i, 'L2')}
)
total_subsequent_iter_error = max(subsequent_iter_error.values())
# check if error criterion has been met
if total_subsequent_iter_error < subdomain.mesh.hmin()/60:
subdomain_calc_isdone[sd_index] = True
# prepare next iteration
# write the newly calculated solution to the inteface dictionaries
# for communication
subdomain.write_pressure_to_interfaces()
for phase in subdomain.has_phases:
subdomain.pressure_prev_iter[phase].vector().set_local(
subdomain.pressure[phase].vector().get_local()
subdomain.pressure_prev_iter[phase].assign(
subdomain.pressure[phase]
)
else:
# one subdomain is done. Check if all subdomains are done to
......@@ -356,6 +363,8 @@ class LDDsimulation(object):
if debug and subdomain.mesh.num_cells() < 36:
print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
## Private methods
def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
mesh_resolution: float = None, #
......@@ -609,7 +618,7 @@ class LDDsimulation(object):
pDCa = df.Expression(pDC[phase], domain = mesh, degree = interpolation_degree, t=time)
self.dirichletBC_dfExpression[num].update({phase: pDCa})
else:
pDCw = self.dirichletBC_dfExpression[num][phase].t = time
pDCa = self.dirichletBC_dfExpression[num][phase].t = time
self.outerBC[num].update(
{phase: df.DirichletBC(V[phase], pDCa, boundary_marker, 1)}
......@@ -636,9 +645,10 @@ class LDDsimulation(object):
f_expression = {
phase: df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
}
# print(f_expression.__dir__())
# subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
subdomain.source = f_expression.copy()
subdomain.source.update(
{phase: f_expression[phase]}#
)
else:
# note that the default interpolation degree is 2
subdomain.source[phase].t = time
......@@ -85,7 +85,7 @@ isRichards = {
}
Tmax = 2
dt1 = 0.01
dt1 = 0.001
dt2 = dt1
# the timestep is taken as a dictionary to be able to run different timesteps on
# different subdomains.
......@@ -185,7 +185,7 @@ write_to_file = {
'L_iterations': True
}
mesh_resolution = 4
mesh_resolution = 100
# initialise LDD simulation class
simulation = ldd.LDDsimulation(tol = 1E-12)
......@@ -239,7 +239,7 @@ df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interfac
# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
simulation.LDDsolver(time = 0)
simulation.LDDsolver(time = 0, debug = True)
# df.info(parameters, True)
# for iterations in range(1):
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment