diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index fb412f3b872ee5114c3cc00d0ee82857fd6fb49c..850cfb5980ed89eb41911256631aa5d91d6af769 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -10,6 +10,12 @@ import mshr
import domainPatch as dp
import helpers as hlp
import pandas as pd
+import sys
+# import errno
+import h5py
+from termcolor import colored
+
+
class LDDsimulation(object):
@@ -147,7 +153,7 @@ class LDDsimulation(object):
'absolute_tolerance': 1E-12,
'relative_tolerance': 1E-10,
'maximum_iterations': 1000,
- 'report': True
+ 'report': False
}
@@ -248,7 +254,32 @@ class LDDsimulation(object):
"Initialising simulation now ... ")
self.initialise()
+ df.info(colored("Start time stepping","green"))
# write initial values to file.
+ while self.t < self.T:
+ print(f"entering timestep t = self.t")
+ # check if the solver should beanalised or not.
+ if analyse_solver_at_times is not None and np.isin(self.t, analyse_solver_at_times):
+ self.LDDsolver(debug = False, analyse_timestep = True)
+ else:
+ self.LDDsolver(debug = False, analyse_timestep = False)
+
+ self.t += self.timestep_size
+ # write solutions to files.
+ for subdom_ind, subdomain in self.subdomain.items():
+ subdomain.write_solution_to_xdmf(#
+ file = self.solution_file[subdom_ind], #
+ time = self.t,#
+ write_iter_for_fixed_time = False,
+ )
+
+ df.info("Finished")
+ df.info("Elapsed time:"+str(timer.elapsed()[0]))
+ timer.stop()
+ # r_cons.close()
+ # energy.close()
+ df.list_timings(df.TimingClear.keep, [df.TimingType.wall, df.TimingType.system])
+ df.dump_timings_to_xml(self.output_dir+"timings.xml",df.TimingClear.keep)
def LDDsolver(self, time: float = None, #
@@ -378,14 +409,13 @@ class LDDsimulation(object):
# write the newly calculated solution to the inteface dictionaries
# for communication
subdomain.write_pressure_to_interfaces()
+ if analyse_timestep:
+ subdomain.write_solution_to_xdmf(#
+ file = solution_over_iteration_within_timestep[sd_index], #
+ time = time,#
+ write_iter_for_fixed_time = True,
+ )#subsequent_errors = subsequent_errors_over_iteration[sd_index])
for phase in subdomain.has_phases:
- if analyse_timestep:
- subdomain.write_solution_to_xdmf(#
- file = solution_over_iteration_within_timestep[sd_index], #
- time = time,#
- write_iter_for_fixed_time = True,
- )#subsequent_errors = subsequent_errors_over_iteration[sd_index])
-
subdomain.pressure_prev_iter[phase].assign(
subdomain.pressure[phase]
)
@@ -720,7 +750,7 @@ class LDDsimulation(object):
lambda_param = self.lambda_param[subdom_num],#
relative_permeability = self.relative_permeability[subdom_num],#
saturation = self.saturation[subdom_num],#
- timestep_size = self.timestep_size[subdom_num],#
+ timestep_size = self.timestep_size,#
tol = self.tol)
}
)
diff --git a/RR-2-patch-test-case/RR-2-patch-test.py b/RR-2-patch-test-case/RR-2-patch-test.py
index 4a4624ead2a738b6dfd30aa99b5bcd714e04a714..c61dff9c4b907843a6272aa7b9711271548d1675 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -84,15 +84,10 @@ isRichards = {
2: True
}
-Tmax = 2
-dt1 = 0.001
-dt2 = dt1
-# the timestep is taken as a dictionary to be able to run different timesteps on
-# different subdomains.
-timestep_size = {
- 1: dt1,
- 2: dt2
-}
+Tmax = 1
+timestep_size = 0.005
+
+time_interval = [0, Tmax]
viscosity = {#
# subdom_num : viscosity
@@ -114,8 +109,8 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
- 1 : {'wetting' :4},#
- 2 : {'wetting' :4}
+ 1 : {'wetting' :100},#
+ 2 : {'wetting' :100}
}
## relative permeabilty functions on subdomain 1
@@ -171,8 +166,18 @@ exact_solution = {
}
# similary to the outer boundary dictionary, if a patch has no outer boundary
-# None should be written instead of an expression.
-dirichletBC = exact_solution
+# None should be written instead of an expression. This is a bit of a brainfuck:
+# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
+# Since a domain patch can have several disjoint outer boundary parts, the expressions
+# need to get an enumaration index which starts at 0. So dirichletBC[ind][j] is
+# the dictionary of outer dirichlet conditions of subdomain ind and boundary part j.
+# finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting'] return
+# the actual expression needed for the dirichlet condition for both phases if present.
+dirichletBC = {
+#subdomain index: {outer boudary part index: {phase: expression}}
+ 1: { 0: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])'}},
+ 2: { 0: {'wetting': '1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])'}}
+}
# def saturation(pressure, subdomain_index):
# # inverse capillary pressure-saturation-relationship
@@ -185,11 +190,11 @@ write_to_file = {
'L_iterations': True
}
-mesh_resolution = 100
+mesh_resolution = 30
# initialise LDD simulation class
simulation = ldd.LDDsimulation(tol = 1E-12)
-simulation.set_parameters(output_dir = "./",#
+simulation.set_parameters(output_dir = "./output/",#
subdomain_def_points = subdomain_def_points,#
isRichards = isRichards,#
interface_def_points = interface_def_points,#
@@ -202,6 +207,7 @@ simulation.set_parameters(output_dir = "./",#
lambda_param = lambda_param,#
relative_permeability = relative_permeability,#
saturation = sat_pressure_relationship,#
+ time_interval = time_interval,#
timestep_size = timestep_size,#
sources = source_expression,#
initial_conditions = initial_condition,#
@@ -213,9 +219,9 @@ simulation.set_parameters(output_dir = "./",#
simulation.initialise()
# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
# subdomain marker functions
-output_dir = simulation.output_dir
-domain_marker = simulation.domain_marker
-mesh_subdomain = simulation.mesh_subdomain
+# output_dir = simulation.output_dir
+# domain_marker = simulation.domain_marker
+# mesh_subdomain = simulation.mesh_subdomain
# mesh = mesh_subdomain[0]
@@ -229,17 +235,18 @@ mesh_subdomain = simulation.mesh_subdomain
# interface.mark(interface_marker, 1)
# simulation._init_interfaces(interface_def_points, adjacent_subdomains)
-interface = simulation.interface
-interface_marker = simulation.interface_marker
-
-subdoms = simulation.subdomain
+# interface = simulation.interface
+# interface_marker = simulation.interface_marker
+#
+# subdoms = simulation.subdomain
-df.File('./subdomain1_interface_marker.pvd') << simulation.subdomain[1].interface_marker
-df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interface_marker
+# df.File('./subdomain1_interface_marker.pvd') << simulation.subdomain[1].interface_marker
+# df.File('./subdomain2_interface_marker.pvd') << simulation.subdomain[2].interface_marker
# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
-
-simulation.LDDsolver(time = 0, debug = True)
+analyse_timesteps = [0, 100*timestep_size, 1, 2, 3]
+simulation.run(analyse_timesteps)
+# simulation.LDDsolver(time = 0, debug = True, analyse_timestep = True)
# df.info(parameters, True)
# for iterations in range(1):