move BoundaryPart and Interface classes to boundary_and_interface.py

LDDsimulation/LDDsimulation.py

@@ -8,6 +8,7 @@ import numpy as np
 import typing as tp
 import mshr
 import domainPatch as dp
+import boundary_and_interface as bi
 import helpers as hlp
 import pandas as pd
 import sys
@@ -15,7 +16,7 @@ import copy
 # import errno
 import h5py
 from termcolor import colored
-
+from solutionFile import SolutionFile
 
 
 
@@ -27,7 +28,9 @@ class LDDsimulation(object):
     """
 
     def __init__(self, dimension: int = 2,#
-                 tol: float = None):
+                 debug: bool = False,#
+                 tol: float = None,#
+                 LDDsolver_tol: float = 1E-8):
         hlp.print_once("\n ###############################################\n",
                        "############# Begin Simulation ################\n",
                        "###############################################\n")
@@ -162,7 +165,8 @@ class LDDsimulation(object):
             'report': False
         }
 
-        self.debug_solver = True
+        self.debug_solver = debug
+        self.LDDsolver_tol = LDDsolver_tol
         # Dictionaries needed by the LDD Solver
         # dictionary holding bool variables for each subdomain indicating wether
         # minimal error has been achieved, i.e. the calculation can be considered
@@ -372,7 +376,7 @@ class LDDsimulation(object):
         # of subsequent iterations of both phases on the subdomain.
         # max(subsequent_error) is used as a measure if the solver is done or not.
         max_subsequent_error = np.zeros(self._number_of_subdomains, dtype=float)
-        error_criterion = 1E-8#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
+        error_criterion = self.LDDsolver_tol#min(1E-9, subdomain.mesh.hmax()**2*subdomain.timestep_size)
         # in case analyse_timestep is None, solution_over_iteration_within_timestep is None
         solution_over_iteration_within_timestep = self.prepare_LDDsolver(analyse_timestep = analyse_timestep)
         ### actual iteration starts here
@@ -402,6 +406,8 @@ class LDDsimulation(object):
                         subsequent_iter_error.update(
                             {phase: error_calculated}
                             )
+                        if debug:
+                            print(f"time = {time}: subsequent error on subdomain {sd_index} for {phase} phase in iteration {iteration} = ", subsequent_iter_error[phase])
 
                     # calculate the maximum over phase of subsequent errors for the
                     # stopping criterion.
@@ -413,7 +419,7 @@ class LDDsimulation(object):
                         subsequent_error_filename = self.output_dir\
                             +self.output_filename_parameter_part[sd_index]\
                             +"subsequent_iteration_errors" +"_at_time"+\
-                            "{number:.{digits}f}".format(number=time, digits=2) +".csv"
+                            "{number:.{digits}f}".format(number=time, digits=4) +".csv"
                         self.write_subsequent_errors_to_csv(
                             filename = subsequent_error_filename, #
                             subdomain_index = sd_index,
@@ -438,7 +444,7 @@ class LDDsimulation(object):
             # end loop over subdomains
             ##### stopping criterion for the solver.
             # only check if error criterion has been met after at least one iteration.
-            for iteration > 1:
+            if iteration > 1:
                 if debug:
                     # print(f" total_subsequent_error in iter {iteration} = {total_subsequent_iter_error }\n")
                     print(f"the maximum distance between any to points in a cell mesh.hmax() on subdomain{sd_index} is h={subdomain.mesh.hmax()}")
@@ -446,15 +452,17 @@ class LDDsimulation(object):
                     print(f"the error criterion is tol = {error_criterion}")
 
                 total_subsequent_error = np.amax(max_subsequent_error)
+                if debug:
+                    print(f"time = {time}: max subsequent error of both subdomain in iteration {iteration} = ", total_subsequent_error)
                 if total_subsequent_error < error_criterion:
                     # if debug:
-                    print(f"subdomain{sd_index} reached error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
+                    print(f"all phases on all subdomains reached a subsequent error = {total_subsequent_error} < tol = {error_criterion} after {iteration} iterations.")
                     all_subdomains_are_done = True
         # end iteration while loop.
 
     def prepare_LDDsolver(self, time: float = None, #
                   debug: bool = False, #
-                  analyse_timestep: bool = False) -> analysing_solution_file_dict: tp.Dict[ind, SolutionFile]:
+                  analyse_timestep: bool = False) -> tp.Dict[int, tp.Type[SolutionFile]]:
         """ initialise LDD iteration for current time time.
 
         This method executes a few prelininary steps before actually starting
@@ -538,6 +546,7 @@ class LDDsimulation(object):
             # method works properly. The current pressure dictionaries should already
             # be populated with the current solution, sinces the solver writes the pressure
             # to interfaces after each iteration.
+            subdomain.write_pressure_to_interfaces()
             subdomain.write_pressure_to_interfaces(previous_iter=True)
             subdomain.calc_gl0_term()
 
@@ -863,7 +872,7 @@ class LDDsimulation(object):
         self.interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
         self.interface_marker.set_all(0)
         for num, vertices in enumerate(interface_def_points):
-            self.interface.append(dp.Interface(vertices=vertices, #
+            self.interface.append(bi.Interface(vertices=vertices, #
                                                 tol=self.tol,#mesh.hmin()/100,
                                                 internal=True,#
                                                 adjacent_subdomains = adjacent_subdomains[num],#
@@ -1138,16 +1147,3 @@ 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}]")
-
-# https://github.com/geo-fluid-dynamics/phaseflow-fenics
-# adapted from Lucas Ostrowski
-class SolutionFile(df.XDMFFile):
-    """
-    This class extends `df.XDMFFile` with some minor changes for convenience.
-    """
-    def __init__(self, mpicomm, filepath):
-        df.XDMFFile.__init__(self, mpicomm, filepath)
-
-        self.parameters["functions_share_mesh"] = True
-        self.parameters["flush_output"] = True
-        self.path = filepath # Mimic the file path attribute from a `file` returned by `open`