From 0df3b2a5982dfdb070b32b6cac603d14c4072a72 Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Sat, 6 Apr 2019 14:45:22 +0200
Subject: [PATCH] fix call by value/call by reference bugs

---
 LDDsimulation/LDDsimulation.py          | 113 +++++++++++++++---------
 LDDsimulation/domainPatch.py            |  62 +++++++------
 RR-2-patch-test-case/RR-2-patch-test.py |  10 ++-
 3 files changed, 115 insertions(+), 70 deletions(-)

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index ba5dc8f..364ff74 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -42,7 +42,7 @@ class LDDsimulation(object):
         df.parameters["form_compiler"]["cpp_optimize"] = True
         # df.parameters['std_out_all_processes'] = False
 
-        df.set_log_level(df.LogLevel.PROGRESS)
+        df.set_log_level(df.LogLevel.INFO)
         # df.set_log_level(df.LogLevel.INFO)
 
         # CRITICAL  = 50, // errors that may lead to data corruption and suchlike
@@ -120,20 +120,8 @@ class LDDsimulation(object):
         # df.list_krylov_solver_preconditioners()
 
         ### Define the linear solver to be used.
-        solver = 'bicgstab' # biconjugate gradient stabilized method
-        preconditioner = 'ilu' # incomplete LU factorization
-        self.linear_solver = df.KrylovSolver(solver, preconditioner)
-        # we use the previous iteration as an initial guess for the linear solver.
-        solver_param = self.linear_solver.parameters
-        solver_param['nonzero_initial_guess'] = True
-        # solver_param['absolute_tolerance'] = 1E-12
-        # solver_param['relative_tolerance'] = 1E-9
-        solver_param['maximum_iterations'] = 1000
-        solver_param['report'] = True
-        # ## print out set solver parameters
-        # for parameter, value in self.linear_solver.parameters.items():
-        #     print(f"parameter: {parameter} = {value}")
-        df.info(solver_param, True)
+        self.solver = 'bicgstab' # biconjugate gradient stabilized method
+        self.preconditioner = 'hypre_amg' # incomplete LU factorization
 
     def set_parameters(self,#
                        output_dir: str,#
@@ -154,7 +142,8 @@ class LDDsimulation(object):
                        sources: tp.Dict[int, tp.Dict[str, str]],#
                        initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
                        dirichletBC_Expressions: tp.Dict[int, tp.Dict[str, str]],#
-                       exact_solution_case: bool = False#
+                       write2file: tp.Dict[str, bool],#
+                       exact_solution_case: bool = False, #
                        )-> None:
 
         """ set parameters of an instance of class LDDsimulation"""
@@ -176,6 +165,7 @@ class LDDsimulation(object):
         self.initial_conditions = initial_conditions
         self.dirichletBC_Expressions = dirichletBC_Expressions
         self.exact_solution_case = exact_solution_case
+        self.write2file = write2file
         self._parameters_set = True
 
     def initialise(self) -> None:
@@ -204,6 +194,21 @@ class LDDsimulation(object):
         # before the iteration gets started all iteration numbers must be nulled
         # and the subdomain_calc_isdone dictionary must be set to False
         for ind, subdomain in self.subdomain.items():
+            #create solver object for subdomains
+            if time < self.calc_tol:
+                subdomain.linear_solver = df.KrylovSolver(self.solver, self.preconditioner)
+                # we use the previous iteration as an initial guess for the linear solver.
+                solver_param = subdomain.linear_solver.parameters
+                solver_param['nonzero_initial_guess'] = True
+                # solver_param['absolute_tolerance'] = 1E-12
+                # solver_param['relative_tolerance'] = 1E-9
+                solver_param['maximum_iterations'] = 1000
+                solver_param['report'] = True
+                # ## print out set solver parameters
+                # for parameter, value in self.linear_solver.parameters.items():
+                #     print(f"parameter: {parameter} = {value}")
+                # df.info(solver_param, True)
+
             subdomain_calc_isdone.update({ind: False})
             # reset all interface[has_interface].current_iteration[ind] = 0, for
             # index has_interface in subdomain.has_interface
@@ -342,15 +347,19 @@ class LDDsimulation(object):
             if debug and subdomain.mesh.num_cells() < 36:
                 print("\npressure before solver:\n", subdomain.pressure[phase].vector().get_local())
 
-            solver = self.linear_solver
-            solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
+            linear_solver = subdomain.linear_solver
+            if linear_solver is None:
+                df.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled, self.solver, self.preconditioner)
+            else:
+                linear_solver.solve(form_assembled, subdomain.pressure[phase].vector(), rhs_assembled)
 
             if debug and subdomain.mesh.num_cells() < 36:
                 print("\npressure after solver:\n", subdomain.pressure[phase].vector().get_local())
 
     ## Private methods
     def _init_meshes_and_markers(self, subdomain_def_points: tp.List[tp.List[df.Point]] = None,#
-                                mesh_resolution: float = None) -> None:
+                                mesh_resolution: float = None, #
+                                ) -> None:
         """ Generate meshes and markerfunctions for sudomains.
 
         This method generate meshes and markerfunctions for sudomains and
@@ -388,6 +397,7 @@ class LDDsimulation(object):
 
         if not mesh_resolution:
             mesh_resolution = self.mesh_resolution
+
         ### generate global mesh and subdomain meshes. Also count number of subdomains.
         # define the global simulation domain polygon first.
         domain = mshr.Polygon(subdomain_def_points[0])
@@ -399,6 +409,11 @@ class LDDsimulation(object):
 
         # create global mesh. Submeshes still need to be extracted.
         mesh = mshr.generate_mesh(domain, mesh_resolution)
+        if self.write2file['meshes_and_markers']:
+            filepath = self.output_dir+"global_mesh.pvd"
+            df.File(filepath) << mesh
+            filepath = self.output_dir+"global_mesh.xml.gz"
+            df.File(filepath) << mesh
         self.mesh_subdomain.append(mesh)
         # define domainmarker on global mesh and set marker values to domain
         # number as assigned in the above for loop.
@@ -460,7 +475,9 @@ class LDDsimulation(object):
                                                 tol=self.tol,#mesh.hmin()/100,
                                                 internal=True,#
                                                 adjacent_subdomains = adjacent_subdomains[num],#
-                                                isRichards = self.isRichards))#
+                                                isRichards = self.isRichards,#
+                                                global_index = num)
+                                  )#
             # the marking number of each interface corresponds to the mathematical
             # numbering of the interfaces (starting with 1 not 0 for the first interface).
             # The reason is that we want to mark outer boundaries with the value 0.
@@ -633,28 +650,42 @@ class LDDsimulation(object):
         """ evaluate time dependent source terms or initialise them if time == 0
         """
         subdomain = self.subdomain[subdomain_index]
-        if np.abs(time) < self.calc_tol:
-            # here t = 0 and we have to initialise the sources.
-            mesh = subdomain.mesh
-            V = subdomain.function_space
-            # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
-            f = self.sources[subdomain_index]
-            if subdomain.isRichards:
+        for phase in subdomain.has_phases:
+            if np.abs(time) < self.calc_tol:
+                # here t = 0 and we have to initialise the sources.
+                mesh = subdomain.mesh
+                # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
+                f = self.sources[subdomain_index][phase]
                 # note that the default interpolation degree is 2
-                fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
+                f_expression = df.Expression(f, domain = mesh, degree = interpolation_degree, t = time)
+                # print(f_expression.__dir__())
                 # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
-                subdomain.source = {'wetting' : fw}
+                subdomain.source.update({phase : f_expression})
             else:
-                fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
-                fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
-                subdomain.source = {'wetting' : fw,#
-                                    'nonwetting' : fnw}
-        else:
-            # here the sources are already initiated and stored in the subdomain.
-            # they must be updated in time.
-            if subdomain.isRichards:
                 # note that the default interpolation degree is 2
-                subdomain.source['wetting'].t = time
-            else:
-                subdomain.source['wetting'].t = time
-                subdomain.source['nonwetting'].t = time
+                subdomain.source[phase].t = time
+
+        # if np.abs(time) < self.calc_tol:
+        #     # here t = 0 and we have to initialise the sources.
+        #     mesh = subdomain.mesh
+        #     # p0 contains both pw_0 and pnw_0 for subdomain[subdomain_index]
+        #     f = self.sources[subdomain_index]
+        #     if subdomain.isRichards:
+        #         # note that the default interpolation degree is 2
+        #         fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
+        #         # subdomain.source = {'wetting' : df.interpolate(fw, V['wetting'])}
+        #         subdomain.source = {'wetting' : fw}
+        #     else:
+        #         fw = df.Expression(f['wetting'], domain = mesh, degree = interpolation_degree, t = time)
+        #         fnw = df.Expression(f['nonwetting'], domain = mesh, degree = interpolation_degree, t = time)
+        #         subdomain.source = {'wetting' : fw,#
+        #                             'nonwetting' : fnw}
+        # else:
+        #     # here the sources are already initiated and stored in the subdomain.
+        #     # they must be updated in time.
+        #     if subdomain.isRichards:
+        #         # note that the default interpolation degree is 2
+        #         subdomain.source['wetting'].t = time
+        #     else:
+        #         subdomain.source['wetting'].t = time
+        #         subdomain.source['nonwetting'].t = time
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 80774e5..32ecc6e 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -161,7 +161,8 @@ class DomainPatch(df.SubDomain):
         self.ds = None
         # solution function(s) for the form set by _init_function_space
         self.pressure = None
-        self.source = None
+        # this is being populated by LDDsimulation._eval_sources()
+        self.source = dict()
         # FEM function holding the previous iteration.
         self.pressure_prev_iter = None
         # FEM function holding the pressures of the previous timestep t_n.
@@ -179,6 +180,9 @@ class DomainPatch(df.SubDomain):
         # dofs of the previous iteration of the neighbouring pressure. This is used
         # to assemble the gli terms in the method self.calc_gli_term().
         self.neighbouring_interface_pressure = None
+        # valiable holding a solver object of type df.KrylovSolver. It is set the
+        # first time LDDsimulation.LDDsolver is run.
+        self.linear_solver = None
 
         # sometimes we need to loop over the phases and the length of the loop is
         # determined by the model we are assuming
@@ -292,7 +296,8 @@ class DomainPatch(df.SubDomain):
                     # assumes a Richards model, we assume constant normalised atmospheric pressure
                     # and no interface term is practically appearing.
                     # interface_forms += (df.Constant(0)*phi_a)*ds(interface)
-                    interface_forms.append((df.Constant(0)*phi_a)*ds(marker_value))
+                    pass
+                    #interface_forms.append((df.Constant(0)*phi_a)*ds(marker_value))
                 else:
                     # interface_forms += (dt*Lambda_a*pa*phi_a)*ds(interface)
                     interface_forms.append((dt*Lambda_a*pa*phi_a)*ds(marker_value))
@@ -355,31 +360,31 @@ class DomainPatch(df.SubDomain):
 
             if iteration == 0:
                 n = df.FacetNormal(self.mesh)
-                pw_prev_iter = self.pressure_prev_timestep
+                p_prev_iter = self.pressure_prev_timestep
                 k = self.relative_permeability
                 S = self.saturation
                 mu = self.viscosity
 
                 if not neighbour_iter_num == 0:
-                    # save the gl term from the neighbour first to use it in the
+                    # copy the gl term from the neighbour first to use it in the
                     # next iteration
                     interface.gli_term_prev[subdomain].update(#
-                        {'wetting': interface.gli_term[neighbour]['wetting']}
+                        {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
                         )
                     # in case we have two-phase and the neighbour is two-phase, two,
                     # we need to Additionally save the nonwetting values.
                     if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
                         interface.gli_term_prev[subdomain].update(#
-                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
                             )
 
                 if self.isRichards:
                     # assuming that we normalised atmospheric pressure to zero,
                     # pc = pw in the Richards case.
-                    pc_prev = pw_prev_iter['wetting']
+                    pc_prev = p_prev_iter['wetting']
                 else:
-                    pw_prev = pw_prev_iter['wetting']
-                    pnw_prev = pw_prev_iter['nonwetting']
+                    pw_prev = p_prev_iter['wetting']
+                    pnw_prev = p_prev_iter['nonwetting']
                     pc_prev = pnw_prev - pw_prev
 
                 for phase in self.has_phases:
@@ -388,7 +393,7 @@ class DomainPatch(df.SubDomain):
                     # relative permeability
                     ka = k[phase]
                     # previous pressure iteration
-                    pa_prev = pw_prev_iter[phase]
+                    pa_prev = p_prev_iter[phase]
                     # testfunction
                     phi_a = self.testfunction[phase]
                     if phase == 'wetting':
@@ -397,7 +402,7 @@ class DomainPatch(df.SubDomain):
                         # TODO: add intrinsic permeabilty!!!!
                         # TODO: add gravitiy term!!!!!
                         flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
-                        gl0 = (df.dot(flux, n) - Lambda[phase]*pa_prev)
+                        gl0 = df.dot(flux, n) - Lambda[phase]*pa_prev
                         gl0_assembled = df.assemble(dt*gl0*phi_a*ds).get_local()
                         # if dofs_in_common_with_other_interfaces[phase].size == 0
                         # there are no dofs that lie on another interface and we can
@@ -424,7 +429,8 @@ class DomainPatch(df.SubDomain):
                             # if the neighbour of our subdomain (with index self.subdomain_index)
                             # assumes a Richards model, we don't have gli terms for the
                             # nonwetting phase and assemble the terms as zero form.
-                            gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
+                            pass
+                            # gli_assembled_tmp[phase] += df.assemble(df.Constant(0)*ds).get_local()
                         else:
                             # in this case the neighbour assumes two-phase flow
                             # and we need to calculte a gl0 torm for the nonwetting
@@ -461,12 +467,12 @@ class DomainPatch(df.SubDomain):
                     # we assume, that the neighbouring patch has done a previous
                     # iteration and we need the value to calculate the gli term with
                     interface.gli_term_prev[subdomain].update(#
-                        {'wetting': interface.gli_term[neighbour]['wetting']}
+                        {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
                         )
 
                     if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
                         interface.gli_term_prev[subdomain].update(#
-                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
                             )
 
                 for phase in self.has_phases:
@@ -519,12 +525,12 @@ class DomainPatch(df.SubDomain):
                     # gli_term_prev has been used by the above and can now safely
                     # be overwritten with the gli_term of the neighbour.
                     interface.gli_term_prev[subdomain].update(#
-                        {'wetting': interface.gli_term[neighbour]['wetting']}
+                        {'wetting': interface.gli_term[neighbour]['wetting'].copy()}
                         )
 
                     if not interface.neighbour_isRichards[subdomain] and not self.isRichards:
                         interface.gli_term_prev[subdomain].update(#
-                            {'nonwetting': interface.gli_term[neighbour]['nonwetting']}
+                            {'nonwetting': interface.gli_term[neighbour]['nonwetting'].copy()}
                             )
                 else:
                     print(f"neighbour_iter_num = {neighbour_iter_num} for neighbour = {neighbour}\n",
@@ -542,7 +548,7 @@ class DomainPatch(df.SubDomain):
         # dictionaries.
         for phase in self.has_phases:
             self.gli_assembled.update(
-                {phase: gli_assembled_tmp[phase]}
+                {phase: gli_assembled_tmp[phase].copy()}
                 )
             for ind in self.has_interface:
                 # self._calc_gli_term_on(interface, iteration)
@@ -682,9 +688,9 @@ class DomainPatch(df.SubDomain):
             )
             for phase in self.has_phases:
                 self._dof_indices_of_interface[ind].update(#
-                    {phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value)}
+                    {phase: self.interface[ind].dofs_on_interface(V[phase], marker, marker_value).copy()}
                 )
-                print(f"subdomain{self.subdomain_index}: dofs on interface {ind}:", self._dof_indices_of_interface[ind][phase])
+                print(f"subdomain{self.subdomain_index}: dofs on interface{ind}:", self._dof_indices_of_interface[ind][phase])
 
     def _calc_interface_has_common_dof_indices(self):
         """ populate dictionary self._interface_has_common_dof_indices
@@ -715,12 +721,12 @@ class DomainPatch(df.SubDomain):
             self._interface_has_common_dof_indices.update({interf_ind: dict()})
             if self.isRichards:
                 self._interface_has_common_dof_indices[interf_ind].update(
-                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]}
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w].copy()}
                 )
             else:
                 self._interface_has_common_dof_indices[interf_ind].update(
-                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w]},
-                    {'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw]}
+                    {'wetting': interf_dofs['wetting'][is_part_of_neighbour_interface_w].copy()},
+                    {'nonwetting': interf_dofs['nonwetting'][is_part_of_neighbour_interface_nw].copy()}
                 )
 
 
@@ -933,7 +939,7 @@ class Interface(BoundaryPart):
                             of whitch this interface is an interface to.
     """
     def __init__(self, #
-                domain_index: int = None,#
+                global_index: int,#
                 adjacent_subdomains: np.array = None,#
                 # this is an array of bools
                 isRichards: np.array = None,#
@@ -941,6 +947,8 @@ class Interface(BoundaryPart):
         #
         self.isRichards = isRichards
         self.adjacent_subdomains = adjacent_subdomains
+        # global index ob the subdomain
+        self.global_index = global_index
         # pressure_values is a dictionary of dictionaries, one for each of the
         # adjacent subdomains, holding the dof values over the interface of the
         # pressure (solution) values for communication. It is initialised by the
@@ -1096,13 +1104,13 @@ class Interface(BoundaryPart):
                 # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
                 # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
                 self.pressure_values[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
-            print("have written pressure interface values\n", self.pressure_values[subdomain_ind][phase])
+            print(f"have written subdomain{subdomain_ind} pressure to interface{self.global_index} for phase {phase}\n", self.pressure_values[subdomain_ind][phase])
         else:
             for dof_index in interface_dofs:
                 # from_function.vector() is orderd by dof and not by vertex numbering of the mesh.
                 # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
                 self.pressure_values_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_function.vector()[dof_index]})
-            print("have written previous pressure interface values\n", self.pressure_values_prev[subdomain_ind][phase])
+            print(f"have written subdomain{subdomain_ind} pressure_prev_iter to interface{self.global_index} for phase {phase}\n", self.pressure_values_prev[subdomain_ind][phase])
 
 
     def read_gli_dofs(self, from_array: np.array, #
@@ -1145,13 +1153,13 @@ class Interface(BoundaryPart):
                 # from_array is orderd by dof and not by vertex numbering of the mesh.
                 # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
                 self.gli_term[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
-            print("have written gl interface values\n", self.gli_term[subdomain_ind][phase])
+            print(f"have written subdomain{subdomain_ind} gl dofs to interface{self.global_index} for phase {phase}\n", self.gli_term[subdomain_ind][phase])
         else:
             for dof_index in interface_dofs:
                 # from_array is orderd by dof and not by vertex numbering of the mesh.
                 # parent needs a mesh vertex index of the submesh, therefore d2v is needed.
                 self.gli_term_prev[subdomain_ind][phase].update({parent[d2v[dof_index]]: from_array[dof_index]})
-            print("have written previous gl interface values\n", self.gli_term_prev[subdomain_ind][phase])
+            print(f"have written subdomain{subdomain_ind} gl_prev dofs to interface{self.global_index} for phase {phase}\n", self.gli_term_prev[subdomain_ind][phase])
 
 
     def write_pressure_dofs(self, to_function: df.Function, #
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 d19ff5e..2b8eeae 100755
--- a/RR-2-patch-test-case/RR-2-patch-test.py
+++ b/RR-2-patch-test-case/RR-2-patch-test.py
@@ -85,7 +85,7 @@ isRichards = {
     }
 
 Tmax = 2
-dt1 = 0.1
+dt1 = 0.01
 dt2 = dt1
 # the timestep is taken as a dictionary to be able to run different timesteps on
 # different subdomains.
@@ -180,11 +180,16 @@ dirichletBC = exact_solution
 #
 # sa
 
+write_to_file = {
+    'meshes_and_markers': True,
+    'L_iterations': True
+}
+
 mesh_resolution = 4
 
 # initialise LDD simulation class
 simulation = ldd.LDDsimulation(tol = 1E-12)
-simulation.set_parameters(output_dir = "",#
+simulation.set_parameters(output_dir = "./",#
     subdomain_def_points = subdomain_def_points,#
     isRichards = isRichards,#
     interface_def_points = interface_def_points,#
@@ -202,6 +207,7 @@ simulation.set_parameters(output_dir = "",#
     initial_conditions = initial_condition,#
     dirichletBC_Expressions = dirichletBC,#
     exact_solution_case = True,#
+    write2file = write_to_file,#
     )
 
 simulation.initialise()
-- 
GitLab