fix meshstudy implementation

LDDsimulation/LDDsimulation.py

@@ -120,8 +120,6 @@ class LDDsimulation(object):
         # self.calc_tol = self.tol
         # list of timesteps that get analyed. Gets initiated by self._init_analyse_timesteps
         self.analyse_timesteps = None
-        # The number of subdomains are counted by self.init_meshes_and_markers()
-        self._number_of_subdomains = 0
         # variable to check if self.set_parameters() has been called
         self._parameters_set = False
         # flag holds true if self.initialise() has been executed.
@@ -257,7 +255,6 @@ class LDDsimulation(object):
         self.saturation = saturation
         # simulation start time and time variable
         self.starttime = starttime
-        self.t = starttime
         # total number of timesteps.
         self.number_of_timesteps = number_of_timesteps
         self.number_of_timesteps_to_analyse = number_of_timesteps_to_analyse
@@ -287,6 +284,8 @@ class LDDsimulation(object):
         """
         if not self._parameters_set:
             raise RuntimeError('Parameters note set!\nYou forgott to run self.set_parameters(**kwds)')
+        # The number of subdomains are counted by self.init_meshes_and_markers()
+        self._number_of_subdomains = 0
         self._add_parameters_to_output_dirname()
         self._init_analyse_timesteps()
         df.info(colored("initialise meshes and marker functions ...\n", "yellow"))
@@ -349,6 +348,7 @@ class LDDsimulation(object):
         df.info(colored("Start time stepping","green"))
         # write initial values to file.
         self.timestep_num = int(1)
+        self.t = self.starttime
 
         # note that at this point of the code self.t = self.starttime as set in
         # the beginning.
@@ -554,7 +554,7 @@ class LDDsimulation(object):
                 # of the subdomain for communication.
                 subdomain.write_pressure_to_interfaces()
 
-                if analyse_timestep:
+                if analyse_timestep and self.write_to_file['L_iterations_per_timestep']:
                     subdomain.write_solution_to_xdmf(#
                         file = solution_over_iteration_within_timestep[sd_index], #
                         time = time,#
@@ -1226,36 +1226,29 @@ class LDDsimulation(object):
                 subdomain.pressure_prev_iter[phase].assign(pa0_interpolated)
                 subdomain.pressure_prev_timestep[phase].assign(pa0_interpolated)
                 # write zeros to the abs_diff variables
+                if self.write2file['absolute_differences']:
+                    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
+                        )
 
-                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
+            if self.write2file['solutions']:
+                subdomain.write_solution_to_xdmf(#
+                    file = self.solution_file[num], #
+                    time = self.starttime,#
+                    write_iter_for_fixed_time = False,
                     )
 
-            # this is also being done in the self.prepare_LDDsolver method!!!!
-            # # populate interface dictionaries with pressure values. The calc_gli_term
-            # # of each subdomain reads from the interface.pressure_values_prev
-            # # and from interface.pressure_values dictionary so both need to be populated.
-            # if self.interface_def_points is not None:
-            #     subdomain.write_pressure_to_interfaces()
-            #     subdomain.write_pressure_to_interfaces(previous_iter = True)
-
-            subdomain.write_solution_to_xdmf(#
-                file = self.solution_file[num], #
-                time = self.starttime,#
-                write_iter_for_fixed_time = False,
-                )
-
 
     def _init_DirichletBC_dictionary(self, interpolation_degree: int = None):
         """initialise dictionary that holds the dolfin.DirichletBC objects
@@ -1494,7 +1487,8 @@ class LDDsimulation(object):
     def _init_simulation_output(self):
         """set up dictionary for simulation output"""
         self.output = dict()
-        for subdom_ind, subdomain in self.subdomain.items():
+        for subdom_ind in self.isRichards.keys():
+            subdomain = self.subdomain[subdom_ind]
             self.output.update({subdom_ind: dict()})
             self.output[subdom_ind].update({'mesh_size': subdomain.mesh_size})
             self.output[subdom_ind].update({'errornorm': dict()})

LDDsimulation/domainPatch.py

@@ -777,7 +777,6 @@ class DomainPatch(df.SubDomain):
         This method sets the class variable
             self.function_space["pressure"]
             self.trialfunction
-            self.pressure
             self.testfunction
 
         INPUT

TP-TP-2-patch-pure-dd-avoid-interface-at-origin/mesh_study_convergence/TP-TP-2-patch-pure-dd-convergence-study.py

@@ -24,16 +24,17 @@ solver_tol = 6E-7
 max_iter_num = 1000
 FEM_Lagrange_degree = 1
 mesh_study = True
+resolutions = [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
 
 ############ GRID #######################
 # mesh_resolution = 20
 timestep_size = 0.0001
-number_of_timesteps = 100
+number_of_timesteps = 200
 plot_timestep_every = 1
 # decide how many timesteps you want analysed. Analysed means, that we write out
 # subsequent errors of the L-iteration within the timestep.
 number_of_timesteps_to_analyse = 0
-starttime = 0
+starttime = 0.0
 
 Lw = 0.25 #/timestep_size
 Lnw=Lw
@@ -47,6 +48,27 @@ analyse_condition = False
 
 output_string = "./output/{}-{}_timesteps{}_P{}-solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)
 
+# toggle what should be written to files
+if mesh_study:
+    write_to_file = {
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': False,
+        'absolute_differences': False,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': False
+    }
+else:
+    write_to_file = {
+        'meshes_and_markers': True,
+        'L_iterations_per_timestep': False,
+        'solutions': True,
+        'absolute_differences': True,
+        'condition_numbers': analyse_condition,
+        'subsequent_errors': True
+    }
+
+
 ##### Domain and Interface ####
 # global simulation domain domain
 sub_domain0_vertices = [df.Point(-1.0,-1.0), #
@@ -513,39 +535,39 @@ for subdomain in isRichards.keys():
 #
 # sa
 
+# toggle what should be written to files
 if mesh_study:
     write_to_file = {
         'meshes_and_markers': True,
-        'L_iterations': True,
+        'L_iterations_per_timestep': False,
         'solutions': False,
         'absolute_differences': False,
-        'condition_numbers': False,
+        'condition_numbers': analyse_condition,
         'subsequent_errors': False
     }
 else:
     write_to_file = {
         'meshes_and_markers': True,
-        'L_iterations': True,
+        'L_iterations_per_timestep': False,
         'solutions': True,
         'absolute_differences': True,
-        'condition_numbers': True,
+        'condition_numbers': analyse_condition,
         'subsequent_errors': True
     }
 
 
-# initialise LDD simulation class
-simulation = ldd.LDDsimulation(
-    tol=1E-14,
-    LDDsolver_tol=solver_tol,
-    debug=debugflag,
-    max_iter_num=max_iter_num,
-    FEM_Lagrange_degree=FEM_Lagrange_degree,
-    mesh_study=mesh_study
-    )
-
-resolutions = [5, 10, 15, 20] #, 20, 25, 30, 35, 40, 45, 50]
-
 for mesh_resolution in resolutions:
+    use_case = use_case + "-mesh-res_{}".format(mesh_resolution)
+    # initialise LDD simulation class
+    simulation = ldd.LDDsimulation(
+        tol=1E-14,
+        LDDsolver_tol=solver_tol,
+        debug=debugflag,
+        max_iter_num=max_iter_num,
+        FEM_Lagrange_degree=FEM_Lagrange_degree,
+        mesh_study=mesh_study
+        )
+
     simulation.set_parameters(use_case=use_case,
                               output_dir=output_string,
                               subdomain_def_points=subdomain_def_points,