add timestepping and write run method

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)
                     }
                 )

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):