add _init_simulation_output method and rewrite output dictionary

LDDsimulation/LDDsimulation.py

@@ -302,6 +302,7 @@ class LDDsimulation(object):
         # solution case.
         self._init_exact_solution_expression()
         self._init_DirichletBC_dictionary()
+        self._init_simulation_output()
         if self.exact_solution and self.write2file['solutions']:
             df.info(colored("writing exact solution for all time steps to xdmf files ...", "yellow"))
             self.write_exact_solution_to_xdmf()
@@ -343,18 +344,7 @@ class LDDsimulation(object):
         df.info(colored("Start time stepping","green"))
         # write initial values to file.
         self.timestep_num = int(1)
-        # set up spacetime_errornorm dictionaries.
-        self.spacetime_errornorm = dict()
-        for subdom_ind, subdomain in self.subdomain.items():
-            self.spacetime_errornorm.update({subdom_ind: dict()})
-            for phase in subdomain.has_phases:
-                self.spacetime_errornorm[subdom_ind].update({phase: dict()})
-                self.spacetime_errornorm[subdom_ind][phase].update(
-                    {# keys mean the norm in time. We use L2 in space
-                    'Linf': 0.0,
-                    'L1': 0.0,
-                    'L2': 0.0 }
-                )
+        
         # note that at this point of the code self.t = self.starttime as set in
         # the beginning.
         while self.timestep_num <= self.number_of_timesteps:
@@ -395,10 +385,10 @@ class LDDsimulation(object):
                         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)
-                        self.spacetime_errornorm[subdom_ind][phase]['L1'] += self.timestep_size*self.spacetime_errornorm[subdom_ind][phase]['L1']
-                        self.spacetime_errornorm[subdom_ind][phase]['L2'] += self.timestep_size*self.spacetime_errornorm[subdom_ind][phase]['L2']**2
-                        if error_calculated > self.spacetime_errornorm[subdom_ind][phase]['Linf']:
-                            self.spacetime_errornorm[subdom_ind][phase]['Linf'] = error_calculated
+                        self.output[subdom_ind]['errornorm'][phase]['L1'] += self.timestep_size*self.output[subdom_ind]['errornorm'][phase]['L1']
+                        self.output[subdom_ind]['errornorm'][phase]['L2'] += self.timestep_size*self.output[subdom_ind]['errornorm'][phase]['L2']**2
+                        if error_calculated > self.output[subdom_ind]['errornorm'][phase]['Linf']:
+                            self.output[subdom_ind]['errornorm'][phase]['Linf'] = error_calculated
 
                         if norm_pa_exact > self.tol:
                             relative_L2_errornorm.update(
@@ -437,7 +427,7 @@ class LDDsimulation(object):
                                 )
 
             self.timestep_num += 1
-        self.spacetime_errornorm[subdom_ind][phase]['L2'] = np.sqrt(self.spacetime_errornorm[subdom_ind][phase]['L2'])
+        self.output[subdom_ind]['errornorm'][phase]['L2'] = np.sqrt(self.output[subdom_ind]['errornorm'][phase]['L2'])
 
         print("LDD simulation use case: {}".format(self.use_case))
         df.info("Finished calculation")
@@ -455,7 +445,7 @@ class LDDsimulation(object):
         # 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"))
-        return self.spacetime_errornorm
+        return self.output
 
     def LDDsolver(self, time: float = None, #
                   debug: bool = False, #
@@ -759,7 +749,7 @@ class LDDsimulation(object):
         for subdom_ind, subdomain in self.subdomain.items():
             tau = subdomain.timestep_size
             L = subdomain.L
-            h = subdomain.mesh.hmax()
+            h = subdomain.mesh_size
             if subdomain.isRichards:
                 Lam = subdomain.lambda_param['wetting']
                 name1 ="subdomain" + "{}".format(subdom_ind)
@@ -1488,3 +1478,20 @@ class LDDsimulation(object):
         self.endtime = self.number_of_timesteps*self.timestep_size
         # time_interval = [starttime, Tmax]
         print(f"The simulation interval is [{self.starttime}, {self.endtime}]")
+        
+        
+    def _init_simulation_output(self):
+        """set up dictionary for simulation output"""
+        self.output = dict()
+        for subdom_ind, subdomain in self.subdomain.items():
+            self.output.update({subdom_ind: dict()})
+            self.output[subdom_ind].update({'mesh_size': subdomain.mesh_size})
+            self.output[subdom_ind].update({'errornorm': dict()})
+            for phase in subdomain.has_phases:
+                self.output[subdom_ind]['errornorm'].update({phase: dict()})
+                self.output[subdom_ind]['errornorm'][phase].update(
+                    {# keys mean the norm in time. We use L2 in space
+                    'Linf': 0.0,
+                    'L1': 0.0,
+                    'L2': 0.0 }
+                )

LDDsimulation/domainPatch.py

@@ -149,6 +149,7 @@ class DomainPatch(df.SubDomain):
         self.subdomain_index = subdomain_index
         self.isRichards = isRichards
         self.mesh = mesh
+        self.mesh_size = self.mesh.hmax()
         self.porosity = porosity
         self.viscosity = viscosity
         self.outer_boundary_def_points = outer_boundary_def_points