fix bugs

LDDsimulation/LDDsimulation.py

@@ -81,9 +81,9 @@ class LDDsimulation(object):
 
         ## Private variables
         # maximal number of L-iterations that the LDD solver uses.
-        self._max_iter_num = 100
+        self._max_iter_num = 5
         # TODO rewrite this with regard to the mesh sizes
-        self.calc_tol = 10^-6
+        self.calc_tol = self.tol
         # 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
@@ -111,6 +111,29 @@ class LDDsimulation(object):
         # These are interpolated and stored to each respective subdomain by
         # self._init_initial_values()
         self.initial_conditions = None
+        ########## SOLVER DEFINITION AND PARAMETERS ########
+        # print("\nLinear Algebra Backends:")
+        # df.list_linear_algebra_backends()
+        # print("\nLinear Solver Methods:")
+        # df.list_linear_solver_methods()
+        # print("\nPeconditioners for Krylov Solvers:")
+        # 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)
 
     def set_parameters(self,#
                        output_dir: str,#
@@ -181,7 +204,7 @@ 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():
-            subdomain_calc_isdone.update({ind, False})
+            subdomain_calc_isdone.update({ind: False})
             # reset all interface[has_interface].current_iteration[ind] = 0, for
             # index has_interface in subdomain.has_interface
             subdomain.null_all_interface_iteration_numbers()
@@ -194,15 +217,15 @@ class LDDsimulation(object):
             # the following only needs to be done when time is not equal 0 since
             # our initialisation functions already took care of setting what follows
             # for the initial iteration step at time = 0.
-            if time > self.tol:
+            if time > self.calc_tol:
                 # this is the beginning of the solver, so iteration must be set
                 # to 0.
                 subdomain.iteration_number = 0
+                for phase in subdomain.has_phases:
                     # set the solution of the previous timestep as new prev_timestep pressure
-                subdomain.pressure_prev_timestep = subdomain.pressure
-                # TODO recheck if
-                subdomain.pressure_prev_iter = subdomain.pressure
-                # needs to be set also
+                    subdomain.pressure_prev_timestep[phase].vector().set_local(
+                        subdomain.pressure[phase].vector().get_local()
+                    )
 
         ### actual iteration starts here
         # gobal stopping criterion for the iteration.
@@ -220,21 +243,24 @@ class LDDsimulation(object):
                     subdomain.iteration_number = iteration
                     # solve the problem on subdomain
                     self.Lsolver_step(subdomain_index = sd_index,#
-                                      debug = True
-                    )
-                    subsequent_iterations_err = self.calc_iteration_error(
-                        subdomain_index = sd_index,#
-                        error_type = "L2"
+                                      debug = True,
                     )
-                    # check if error criterion has been met
-                    if subsequent_iterations_err < self.calc_tol:
-                        subdomain_calc_isdone[sd_index] = True
+                    # subsequent_iterations_err = self.calc_iteration_error(
+                    #     subdomain_index = sd_index,#
+                    #     error_type = "L2"
+                    # )
+                    # # check if error criterion has been met
+                    # if subsequent_iterations_err < self.calc_tol:
+                    #     subdomain_calc_isdone[sd_index] = True
 
                     # prepare next iteration
                     # write the newly calculated solution to the inteface dictionaries
                     # for communication
                     subdomain.write_pressure_to_interfaces()
-                    subdomain.pressure_prev_iter = subdomain.pressure
+                    for phase in subdomain.has_phases:
+                        subdomain.pressure_prev_iter[phase].vector().set_local(
+                            subdomain.pressure[phase].vector().get_local()
+                        )
                 else:
                     # one subdomain is done. Check if all subdomains are done to
                     # stop the while loop if needed.
@@ -255,7 +281,7 @@ class LDDsimulation(object):
             iteration += 1
             # end iteration while loop.
 
-    def Lsolver_step(self, subdomain_index: int],#
+    def Lsolver_step(self, subdomain_index: int,#
                      debug: bool = False) -> None:
         """ L-scheme solver iteration step for an object of class subdomain
 
@@ -273,28 +299,54 @@ class LDDsimulation(object):
         """
         subdomain = self.subdomain[subdomain_index]
         iteration = subdomain.iteration_number
-        if subdomain.isRichards:
+        # calculate the gli terms on the right hand side  and store
+        subdomain.calc_gli_term()
+
+        for phase in subdomain.has_phases:
             # extract L-scheme form and rhs (without gli term) from subdomain.
-            governing_problem = subdomain.governing_problem(phase = 'wetting')
-            form = governing_problem[form]
-            rhs_without_gli = governing_problem[rhs_without_gli]
+            governing_problem = subdomain.governing_problem(phase = phase)
+            form = governing_problem['form']
+            rhs_without_gli = governing_problem['rhs_without_gli']
             # assemble the form and rhs
             form_assembled = df.assemble(form)
             rhs_without_gli_assembled = df.assemble(rhs_without_gli)
-            # calculate the gli term
-            gli_term_assembled = subdomain.calc_gli_term(iteration = iteration)
+            # access the assembled gli term for the rhs.
+            gli_term_assembled = subdomain.gli_assembled[phase]
             # subdomain.calc_gli_term() asslembles gli but on the left hand side
             # so gli_term_assembled needs to actually be subtracted from the rhs.
-            rhs_assembled = rhs_without_gli_assembled - gli_term_assembled
+            if debug and subdomain.mesh.num_cells() < 36:
+                print("\nSystem before applying outer boundary conditions:")
+                print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+            rhs_without_gli_assembled.set_local(rhs_without_gli_assembled.get_local() - gli_term_assembled)
+            rhs_assembled = rhs_without_gli_assembled
+
+            if debug and subdomain.mesh.num_cells() < 36:
+                # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+                # print(f"phase = {phase}: rhs_without_gli:\n", rhs_without_gli_assembled.get_local())
+                print(f"phase = {phase}: gli_term:\n", gli_term_assembled)
+                print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
             # apply outer Dirichlet boundary conditions if present.
             if subdomain.outer_boundary is not None:
-                self.outerBC[subdomain_index]['wetting'].apply(form_assembled, rhs_assembled)
-            solution_pw = subdomain.pressure['wetting'].vector()
+                self.outerBC[subdomain_index][phase].apply(form_assembled, rhs_assembled)
 
-        else:
-            print("implement two phase assembly")
+            if debug and subdomain.mesh.num_cells() < 36:
+                print("\nSystem after applying outer boundary conditions:")
+                # print(f"phase = {phase}: form_assembled:\n", form_assembled.array())
+                print(f"phase = {phase}: rhs_assembled:\n", rhs_assembled.get_local())
+
+            # # set previous iteration as initial guess for the linear solver
+            # pi_prev = subdomain.pressure_prev_iter[phase].vector().get_local()
+            # subdomain.pressure[phase].vector().set_local(pi_prev)
+
+            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)
 
+            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,#
@@ -358,7 +410,7 @@ class LDDsimulation(object):
 
     def _init_interfaces(self, interface_def_points: tp.List[tp.List[df.Point]] = None,#
                               adjacent_subdomains: tp.List[np.ndarray] = None) -> None:
-        """ generate interfaces and interface markerfunctions for sudomains and
+        """ generate interfaces and interface markerfunctions for suddomains and
             set the internal lists used by the LDDsimulation class.
 
         This initialises a list of interfaces and global interface marker functions
@@ -405,7 +457,7 @@ class LDDsimulation(object):
         self.interface_marker.set_all(0)
         for num, vertices in enumerate(interface_def_points):
             self.interface.append(dp.Interface(vertices=vertices, #
-                                                tol=mesh.hmin()/100,#
+                                                tol=self.tol,#mesh.hmin()/100,
                                                 internal=True,#
                                                 adjacent_subdomains = adjacent_subdomains[num],#
                                                 isRichards = self.isRichards))#
@@ -423,6 +475,7 @@ class LDDsimulation(object):
         # Therefor it is used in the subdomain initialisation loop.
         for subdom_num, isR in self.isRichards.items():
             interface_list = self._get_interfaces(subdom_num)
+            # print(f"has_interface list for subdomain {subdom_num}:", interface_list)
             self.subdomain.update(#
                     {subdom_num : dp.DomainPatch(#
                         subdomain_index = subdom_num,#
@@ -438,6 +491,7 @@ class LDDsimulation(object):
                         relative_permeability = self.relative_permeability[subdom_num],#
                         saturation = self.saturation[subdom_num],#
                         timestep_size = self.timestep_size[subdom_num],#
+                        tol = self.tol
                     )})
 
 
@@ -450,7 +504,7 @@ class LDDsimulation(object):
         This method determins all (global) indices of interfaces belonging to
         subdomain with index subdomain_index from adjacent_subdomains.
         """
-        if not adjacent_subdomains:
+        if adjacent_subdomains is None:
             adjacent_subdomains = self.adjacent_subdomains
         else:
             # overwrite what has been set by init()
@@ -480,24 +534,45 @@ class LDDsimulation(object):
             if subdomain.isRichards:
                 # note that the default interpolation degree is 2
                 pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
+                subdomain.pressure = {
+                    'wetting': df.Function(V['wetting'])
+                    }
+                subdomain.pressure_prev_timestep = {
+                    'wetting': df.Function(V['wetting'])
+                    }
                 subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting'])}
-                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting'])}
+                # print("vector()",subdomain.pressure_prev_iter['wetting'].vector().get_local())
+                pw_prev = subdomain.pressure_prev_iter['wetting'].vector().get_local()
+                subdomain.pressure_prev_timestep['wetting'].vector().set_local(pw_prev)
+                subdomain.pressure['wetting'].vector().set_local(pw_prev)
 
             else:
                 pw0 = df.Expression(p0['wetting'], domain = mesh, degree = interpolation_degree)
                 pnw0 = df.Expression(p0['nonwetting'], domain = mesh, degree = interpolation_degree)
+                subdomain.pressure = {
+                    'wetting': df.Function(V['wetting']),#
+                    'nonwetting': df.Function(V['nonwetting'])
+                    }
+                subdomain.pressure_prev_timestep = {
+                    'wetting': df.Function(V['wetting']),
+                    'nonwetting': df.Function(V['nonwetting'])
+                    }
                 subdomain.pressure_prev_iter = {'wetting' : df.interpolate(pw0, V['wetting']),#
                                                 'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
-                subdomain.pressure_prev_timestep = {'wetting' : df.interpolate(pw0, V['wetting']),#
-                                                    'nonwetting' : df.interpolate(pnw0, V['nonwetting'])}
+
+                pw_prev = subdomain.pressure_prev_iter['wetting'].vector().get_local()
+                pnw_prev = subdomain.pressure_prev_iter['nonwetting'].vector().get_local()
+                subdomain.pressure_prev_timestep['wetting'].vector().set_local(pw_prev)
+                subdomain.pressure['wetting'].vector().set_local(pw_prev)
+                subdomain.pressure_prev_timestep['nonwetting'].vector().set_local(pnw_prev)
+                subdomain.pressure['nonwetting'].vector().set_local(pnw_prev)
             # populate interface dictionaries with pressure values.
             subdomain.write_pressure_to_interfaces(initial_values = True)
 
-    def update_DirichletBC_dictionary(self,
+    def update_DirichletBC_dictionary(self, subdomain_index: int,#
                         interpolation_degree: int = 2, #
                         time: float = 0,#
-                        # exact_solution: bool = False,#
-                        subdomain_index: int):
+                        ):
         """ update time of the dirichlet boundary object for subdomain with
             index subdomain_index.
 
@@ -512,7 +587,7 @@ class LDDsimulation(object):
             mesh = subdomain.mesh
             V = subdomain.function_space
             # Here the dictionary has to be created first because t = 0.
-            if np.abs(time) < self.tol :
+            if np.abs(time) < self.calc_tol:
                 self.outerBC.update({num: dict()})
                 self.dirichletBC_dfExpression.update({num: dict()})
 
@@ -529,7 +604,7 @@ class LDDsimulation(object):
 
                 self.outerBC[num].update({'wetting': df.DirichletBC(V['wetting'], pDCw, boundary_marker, 1)})
             else:
-                if np.abs(time) < self.tol :
+                if np.abs(time) < self.calc_tol :
                     # time = 0 so the Dolfin Expression has to be created first.
                     pDCw = df.Expression(pDC['wetting'], domain = mesh, degree = interpolation_degree, t=time)
                     pDCnw = df.Expression(pDC['nonwetting'], domain = mesh, degree = interpolation_degree, t=time)
@@ -550,13 +625,14 @@ class LDDsimulation(object):
             print("subdomain is an inner subdomain and has no outer boundary.\n",
             "no Dirichlet boundary has been set.")
 
-    def _eval_sources(self, interpolation_degree: int = 2, #
-                            time: float = 0,
-                            subdomain_index: int):
+    def _eval_sources(self, subdomain_index: int, #
+                            interpolation_degree: int = 2, #
+                            time: float = 0,#
+                            ):
         """ evaluate time dependent source terms or initialise them if time == 0
         """
         subdomain = self.subdomain[subdomain_index]
-        if np.abs(time) < self.tol:
+        if np.abs(time) < self.calc_tol:
             # here t = 0 and we have to initialise the sources.
             mesh = subdomain.mesh
             V = subdomain.function_space

RR-2-patch-test-case/RR-2-patch-test.py

@@ -183,7 +183,7 @@ dirichletBC = exact_solution
 mesh_resolution = 3
 
 # initialise LDD simulation class
-simulation = ldd.LDDsimulation()
+simulation = ldd.LDDsimulation(tol = 1E-12)
 simulation.set_parameters(output_dir = "",#
     subdomain_def_points = subdomain_def_points,#
     isRichards = isRichards,#
@@ -211,7 +211,7 @@ domain_marker = simulation.domain_marker
 mesh_subdomain = simulation.mesh_subdomain
 
 # mesh = mesh_subdomain[0]
-# interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
+# interface_marker = df.MeshFunction('int', mesh, mesh.topology().dim()-1)
 # interface_marker.set_all(0)
 # interface = dp.Interface(vertices=interface_def_points[0], #
 #                                     tol=mesh.hmin()/100,#
@@ -226,53 +226,68 @@ interface_marker = simulation.interface_marker
 
 subdoms = simulation.subdomain
 
-for iterations in range(1):
-    #
-
-    a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
-    rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
-    rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
-    extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
-    splitrhs = df.assemble(rhssplit1)
-    print("splitrhs: \n", splitrhs.get_local())
-    extrat_assembled = df.assemble(extratrem)
-    print("extratrem: \n", extrat_assembled.get_local())
-    added = splitrhs + extrat_assembled
-    print("both added together:\n", added.get_local())
-
-    A1 = df.assemble(a1)
-    dsform = (timestep*lambda1*u1*v1)*ds1(1)
-    dsform_vec = df.assemble(dsform)
-    b1 = df.assemble(rhs1)
-    #p_gamma1.apply(A1,b1)
-    # outerBC1.apply(A1,b1)
-    u  = df.Function(V1)
-    U1 = u.vector()
-    df.solve(A1,U1,b1)
-    # print("Form:\n", A1.array())
-    print("RHS:\n", b1.get_local())
-    #print("ds term:\n", dsform_vec.array())
-    # print('solution:\n', U1.get_local())
-    p1_i.vector()[:] = U1
-    interface[0].read_pressure_dofs(from_function = p1_i, #
-                            interface_dofs = dofs_on_interface1,#
-                            dof_to_vert_map = dom1_d2v,#
-                            local_to_parent_vertex_map = dom1_parent_vertex_indices,#
-                            phase = 'wetting',#
-                            subdomain_ind = 1)
-#
-# Save mesh to file
-df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
-df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
-df.File('./test_global_interface_marker.pvd') << interface_marker
-df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
-df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
-df.File('./test_domain_marker.pvd') << domain_marker
-df.File('./test_domain_layered_soil_solution.pvd') << u
-df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
-df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
-df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
-df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
+# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
+df.File('./global_interface_marker.pvd') << interface_marker
+# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
+# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
+df.File('./domain_marker.pvd') << domain_marker
+# df.File('./test_domain_layered_soil_solution.pvd') << u
+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)
+# df.info(parameters, True)
+
+# for iterations in range(1):
+#     #
+#
+#     a1 = (L1*u1*v1)*dx1 + (df.inner(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx_1 + (timestep*lambda1*u1*v1)*ds1(1) #timestep*
+#     rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1(1)
+#     rhssplit1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1
+#     extratrem = (timestep*(source1 - g1_i)*v1)*ds1(1)
+#     splitrhs = df.assemble(rhssplit1)
+#     print("splitrhs: \n", splitrhs.get_local())
+#     extrat_assembled = df.assemble(extratrem)
+#     print("extratrem: \n", extrat_assembled.get_local())
+#     added = splitrhs + extrat_assembled
+#     print("both added together:\n", added.get_local())
+#
+#     A1 = df.assemble(a1)
+#     dsform = (timestep*lambda1*u1*v1)*ds1(1)
+#     dsform_vec = df.assemble(dsform)
+#     b1 = df.assemble(rhs1)
+#     #p_gamma1.apply(A1,b1)
+#     # outerBC1.apply(A1,b1)
+#     u  = df.Function(V1)
+#     U1 = u.vector()
+#     df.solve(A1,U1,b1)
+#     # print("Form:\n", A1.array())
+#     print("RHS:\n", b1.get_local())
+#     #print("ds term:\n", dsform_vec.array())
+#     # print('solution:\n', U1.get_local())
+#     p1_i.vector()[:] = U1
+#     interface[0].read_pressure_dofs(from_function = p1_i, #
+#                             interface_dofs = dofs_on_interface1,#
+#                             dof_to_vert_map = dom1_d2v,#
+#                             local_to_parent_vertex_map = dom1_parent_vertex_indices,#
+#                             phase = 'wetting',#
+#                             subdomain_ind = 1)
+# #
+# # Save mesh to file
+# df.File('./test_domain_layered_soil.xml.gz') << mesh_subdomain[0]
+# df.File('./test_domain_mesh.pvd') << mesh_subdomain[0]
+# df.File('./test_global_interface_marker.pvd') << interface_marker
+# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
+# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
+# df.File('./test_domain_marker.pvd') << domain_marker
+# df.File('./test_domain_layered_soil_solution.pvd') << u
+# df.File('./test_subdomain1_interface_marker.pvd') << interface_marker1
+# df.File('./test_subdomain2_interface_marker.pvd') << interface_marker2
+# df.File('./test_subdomain1_boundary_marker.pvd') << boundary_marker1
+# df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
 
 # boundary_marker1 = simulation.subdomain[1].outer_boundary_marker
 # boundary_marker2 = simulation.subdomain[2].outer_boundary_marker
@@ -281,8 +296,8 @@ df.File('./test_subdomain2_boundary_marker.pvd') << boundary_marker2
 # ds_1 = simulation.subdomain[1].ds
 # ds_2 = simulation.subdomain[2].ds
 #
-# interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
-# interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
+# interface_marker1 = simulation.subdomain[1].interface_marker#df.MeshFunction('int', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
+# interface_marker2 = simulation.subdomain[2].interface_marker#df.MeshFunction('int', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
 # # interface_marker1.set_all(0)
 # # interface_marker2.set_all(0)
 # # print(dir(interface_marker1))