From 5e4dc21246fb00f4c291c9de13beaf99dd9aab6a Mon Sep 17 00:00:00 2001
From: David Seus <david.seus@ians.uni-stuttgart.de>
Date: Sat, 26 Oct 2019 16:50:04 +0200
Subject: [PATCH] revert implementing extra gravity flux term

---
 LDDsimulation/LDDsimulation.py                |  23 +-
 LDDsimulation/boundary_and_interface.py       |   9 +-
 LDDsimulation/domainPatch.py                  | 280 +++---------------
 .../TP-R-2-patch-pure-dd-realistic.py         |  45 +--
 .../TP-R-2-patch-pure-dd.py                   |  28 +-
 .../injection/TP-TP-2-patch-injection.py      |  38 ++-
 6 files changed, 101 insertions(+), 322 deletions(-)
 mode change 100644 => 100755 Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py

diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 4699e0e..41c91bb 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -382,7 +382,10 @@ class LDDsimulation(object):
                                analyse_condition=False)
 
             # calculate spacetime errornorms for self.output.
-            space_errornorm = self._spacetime_errornorm_calc_step()
+            if self.exact_solution is not None:
+                space_errornorm = self._spacetime_errornorm_calc_step()
+            else:
+                space_errornorm = None
 
             # write solutions to files. It is paramount, that self.timestep_num
             # be only updated after calling self.write_data_to_disc()
@@ -535,13 +538,6 @@ class LDDsimulation(object):
                     subdomain.pressure_prev_iter[phase].assign(
                         subdomain.pressure[phase]
                     )
-                # onyl becomes active in TPR case
-                if self.include_gravity and subdomain.isRichards and subdomain.TPR_occures:
-                        for interface in subdomain.has_interface_with_TP_neighbour:
-                            subdomain.calc_gravity_neumann_flux(
-                                interface_index=interface,
-                                phase="nonwetting"
-                                )
             # end loop over subdomains
             ##### stopping criterion for the solver.
             total_subsequent_error = np.amax(max_subsequent_error)
@@ -652,17 +648,6 @@ class LDDsimulation(object):
                 subdomain.write_pressure_to_interfaces()
                 subdomain.write_pressure_to_interfaces(previous_iter=True)
                 subdomain.calc_gl0_term()
-                # If gravity is included and a Richards subdomain has a TPR coupling,
-                # we need to calculate the Neumann flux of the nonwetting phase
-                # that appears due to gravity
-                if self.include_gravity:
-                    if subdomain.isRichards and subdomain.TPR_occures:
-                        for interface in subdomain.has_interface_with_TP_neighbour:
-                            subdomain.calc_gravity_neumann_flux(
-                                initial=True,
-                                interface_index=interface,
-                                phase="nonwetting"
-                                )
 
         return solution_over_iteration_within_timestep
 
diff --git a/LDDsimulation/boundary_and_interface.py b/LDDsimulation/boundary_and_interface.py
index 7d1a962..121515b 100644
--- a/LDDsimulation/boundary_and_interface.py
+++ b/LDDsimulation/boundary_and_interface.py
@@ -615,8 +615,6 @@ class Interface(BoundaryPart):
             else:  #function == "gli":
                 self.gli_term.update({subdomain_index: dict()})
                 self.gli_term_prev.update({subdomain_index: dict()})
-                self.gravity_neumann_flux.update({subdomain_index: dict()})
-                self.gravity_neumann_flux_prev.update({subdomain_index: dict()})
 
             for phase in has_phases:
                 if function == "pressure":
@@ -625,8 +623,7 @@ class Interface(BoundaryPart):
                 else:  #function == "gli":
                     self.gli_term[subdomain_index].update({phase: dict()})
                     self.gli_term_prev[subdomain_index].update({phase: dict()})
-                    self.gravity_neumann_flux[subdomain_index].update({phase: dict()})
-                    self.gravity_neumann_flux_prev[subdomain_index].update({phase: dict()})
+
                 space = function_space[function][phase]
                 dof_coordinates = self.dofs_and_coordinates(space, interface_marker, interface_marker_value)
                 # the dof dictionaries are grouped by the global (with respect to the
@@ -640,8 +637,6 @@ class Interface(BoundaryPart):
                     else:  #function == "gli":
                         self.gli_term[subdomain_index][phase].update({global_mesh_facet_index: dict()})
                         self.gli_term_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
-                        self.gravity_neumann_flux[subdomain_index][phase].update({global_mesh_facet_index: dict()})
-                        self.gravity_neumann_flux_prev[subdomain_index][phase].update({global_mesh_facet_index: dict()})
 
                     for dof_coord in facet_dof_coord_dict.values():
                         # we want to use the coordinates as key in a dictionary.
@@ -655,8 +650,6 @@ class Interface(BoundaryPart):
                         else:  #function == "gli":
                             self.gli_term[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                             self.gli_term_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
-                            self.gravity_neumann_flux[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
-                            self.gravity_neumann_flux_prev[subdomain_index][phase][global_mesh_facet_index].update({coord_tupel: 0.0})
                     # end init_interface_values
 
 
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index da2596a..248b432 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -248,14 +248,6 @@ class DomainPatch(df.SubDomain):
             self._calc_interface_dof_indices_and_coordinates()
             self._calc_corresponding_dof_indices()
             self._init_interface_values()
-            # # Dictionary of FEM function holding the same dofs gli_assembled above
-            # # but as a df.Function. This is needed for writing out the gli terms when
-            # # analyising the convergence of a timestep.
-            # self.gli_function = dict()
-            # for phase in self.has_phases:
-            #     self.gli_function.update(
-            #         {phase: df.Function(self.function_space["gli"][phase])}
-            #     )
 
     # END constructor
 
@@ -373,16 +365,6 @@ class DomainPatch(df.SubDomain):
             # gather interface forms
             interface_forms = []
             gli_forms = []
-            # the following needs to be executed on a TP domain with R neighbours
-            # in case gravity is included.
-            if self.include_gravity and phase=="nonwetting" and self.TPR_occures:
-                rhs_gravity_term_interface = []
-                for interface in self.has_interface_with_Richards_neighbour:
-                    marker_value = self.interface[interface].marker_value
-                    gravity_neumann_flux = self.read_gravity_neumann_flux(interface_index=interface, phase="nonwetting")
-                    rhs_gravity_term_interface.append(
-                        -dt*gravity_neumann_flux*phi_a*ds(marker_value)
-                    )
             # the marker values of the interfacese are global interfaces index + 1
             for interface in self.has_interface:
                 # print(f"subdomain{self.subdomain_index} has interfaces: {interface}\n")
@@ -427,12 +409,8 @@ class DomainPatch(df.SubDomain):
         if self.has_interface is not None:
             form = form1 + form2 + sum(interface_forms)
             if self.include_gravity:
-                if phase=="nonwetting" and self.TPR_occures:
-                    rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) \
-                    + rhs_gravity + sum(rhs_gravity_term_interface)
-                else:
-                    # z_a included the minus sign already, see self._calc_gravity_expressions
-                    rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
+                # z_a included the minus sign already, see self._calc_gravity_expressions
+                rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms) + rhs_gravity
             else:
                 rhs = rhs1 + rhs2 + dt*source_a*phi_a*dx - sum(gli_forms)
         else:
@@ -449,198 +427,6 @@ class DomainPatch(df.SubDomain):
         return form_and_rhs
 
 
-    def calc_gravity_neumann_flux(self,
-                      interface_index: int,
-                      phase: str,
-                      initial: bool = False,
-                      debug: bool = False
-                      ) -> None:
-        """calculate the gravity term of the neumann flux
-
-        This method is written such that it can calculate the gravity neumann
-        flux for both phases in case it is needed, but is actually only used
-        for the flux of the nonwetting gravity term of a Richards domain with
-        TPR coupling.
-        """
-        subdomain = self.subdomain_index
-        dt = self.timestep_size
-        Lambda = self.lambda_param
-        interf_ind = interface_index
-
-        interface = self.interface[interf_ind]
-
-        if initial:
-            pressure_previous = self.pressure_prev_timestep
-        else:
-            pressure_previous = self.pressure_prev_iter
-
-        if self.isRichards:
-            # assuming that we normalised atmospheric pressure to zero,
-            # pc = -pw in the Richards case.
-            pc_prev_iter = -pressure_previous['wetting']
-        else:
-            pw_prev_iter = pressure_previous['wetting']
-            pnw_prev_iter = pressure_previous['nonwetting']
-            pc_prev_iter = pnw_prev - pw_prev
-
-        # n = df.FacetNormal(self.mesh)
-        S = self.saturation
-
-        # viscosity
-        mu_a = self.viscosity[phase]
-        # relative permeability
-        ka = self.relative_permeability[phase]
-        z_a = self.gravity_term[phase]
-
-        if phase == 'wetting':
-            # z_a included the minus sign, see self._calc_gravity_expressions
-            flux = 1/mu_a*ka(S(pc_prev_iter))*df.grad(z_a)
-        else:
-            # phase == 'nonwetting'
-            flux = 1/mu_a*ka(1-S(pc_prev_iter))*df.grad(z_a)
-
-        flux_function = df.project(flux, self.function_space["flux"][phase])
-        flux_x, flux_y = flux_function.split()
-        gravity_neumann_flux = df.Function(self.function_space["gli"][phase])
-        # # get dictionaries of global dof indices and coordinates along
-        # # the interface. These dictionaries have the facet indices of
-        # # facets belonging to the interface as keys (with respect to
-        # # the submesh numbering) and the dictionary
-        # # containing pairs {global_dof_index: dof_coordinate} for all
-        # # dofs along the facet as values.
-        neumann_flux_dof_coords = self.interface_dof_indices_and_coordinates["gli"][interf_ind][phase]
-        # p_dof_coordinates = self.interface_dof_indices_and_coordinates["pressure"][interf_ind][phase]
-        # flux_dof_coordinates = self.interface_dof_indices_and_coordinates["flux"][interf_ind][phase]
-        # # the attributes local and global for facet numbering mean
-        # # the following: local facet index is the index of the facet
-        # # with respect to the numbering of the submesh belonging to
-        # # the subdomain. global facet index refers to the facet numbering
-        # # of the mesh of the whole computational domain.
-        local2global_facet = interface.local_to_global_facet_map[subdomain]
-        corresponding_dof_index = self.interface_corresponding_dof_index[interf_ind][phase]
-        for local_facet_ind, normal in interface.outer_normals[subdomain].items():
-            neumann_flux_dofs_along_facet = neumann_flux_dof_coords[local_facet_ind]
-            global_facet_ind = local2global_facet[local_facet_ind]
-            for neumann_flux_dof_index, neumann_flux_dof_coord in neumann_flux_dofs_along_facet.items():
-                flux_x_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_x"]
-                flux_y_dof_index = corresponding_dof_index[local_facet_ind][neumann_flux_dof_index]["flux_y"]
-
-                gravity_neumann_flux.vector()[neumann_flux_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
-                    + flux_y.vector()[flux_y_dof_index]*normal[1]
-
-                # save this newly calculated dof to the interface
-                # dictionary neumann_flux_save_dict
-                dof_coord_tupel = (neumann_flux_dof_coord[0], neumann_flux_dof_coord[1])
-                if initial:
-                    # because we don't know if the Richards domain or the TP domain
-                    # will be up next in the solver, we need to save the initial
-                    # gravity neumann flux to both the current and the prev dictionary.
-                    interface.gravity_neumann_flux_prev[subdomain][phase][global_facet_ind].update(
-                        {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
-                        )
-                    interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
-                        {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
-                        )
-                else:
-                    # in the noninitial case we always save to the current dictionary
-                    # similarly to the pressures.
-                    interface.gravity_neumann_flux[subdomain][phase][global_facet_ind].update(
-                        {dof_coord_tupel: gravity_neumann_flux.vector()[neumann_flux_dof_index]}
-                        )
-
-
-    def read_gravity_neumann_flux(self,
-                      interface_index: int,
-                      phase: str,
-                      debug: bool = False) -> df.Function: # tp.Dict[str, fl.Form]
-        """in the case of a TPR coupling read the additional gravity neumann flux
-        for the nonwetting phase from the neighbour
-
-        IMPLEMENTATION
-        same logic as in the calc_gli_term method applies.
-
-        PARAMETERS
-        ------------------------------------
-        interface_index    int       interface index for which to calculate the
-                                     gravity_neumann_flux term.
-        phase              str       phase to calculate the gravity_neumann_flux term for.
-        # DEBUG:           bool      toggle the debug output.
-        """
-        # this is the number of the current iteration of the subdomain.
-        iteration = self.iteration_number
-        subdomain = self.subdomain_index
-        interface = self.interface[interface_index]
-        # gravity_neumann_flux should be initialised by zero.
-        gravity_neumann_flux = df.Function(self.function_space["gli"][phase])
-
-        # neighbour index
-        neighbour = interface.neighbour[subdomain]
-        # update the current_iteration number of subdomain stored in the interface
-        # to the current iteration number of the subdomain.
-        interface.current_iteration[subdomain] = iteration
-        # save the iteration number of the neighbour
-        neighbour_iter_num = interface.current_iteration[neighbour]
-
-        # read previous neighbouring gravity neumann flux from interface.
-        # depending on our iteration number and the iteration number of
-        # the neighbour, weed to either read from the previous pressure
-        # dictionary or the current one.
-        if iteration == neighbour_iter_num + 1:
-            # in this case the neighbour has not yet iterated beyond
-            # the iteration we are currently in, so
-            # interface.gravity_neumann_flux[neighbour] holds the "previous neumann gravity flux"
-            # we need to read these values from the interface.
-            gravity_flux_read_dict = interface.gravity_neumann_flux[neighbour][phase]
-
-        else:
-            # iteration == neighbour_iter_num:
-            # this is the only possible other case. In this case the
-            # neighbour has iterated once beyond the iteration we are
-            # currently in, so interface.gravity_neumann_flux_prev holds the "previous neumann gravity flux"
-            # we need to read from the interface.
-            gravity_flux_read_dict = interface.gravity_neumann_flux_prev[neighbour][phase]
-
-        flux_interface_dofs = self.interface_dof_indices_and_coordinates["gli"][interface_index][phase]
-        local2global_facet = interface.local_to_global_facet_map[subdomain]
-        for local_facet_index, flux_facet_dofs2coordinates in flux_interface_dofs.items():
-            global_facet_index = local2global_facet[local_facet_index]
-            # read gravity_neumann_flux_prev term from the neighbour.
-            neighbour_flux_on_facet = gravity_flux_read_dict[global_facet_index]
-
-            for flux_dof_index, flux_dof_coord in flux_facet_dofs2coordinates.items():
-                # same with the gravity_neumann_flux_previous_dof
-                gravity_neumann_flux_neighbour_dof = None
-                for flux_neighbour_coord_tupel, flux_neighbour_dof in neighbour_flux_on_facet.items():
-                    flux_neighbour_coord = np.array([flux_neighbour_coord_tupel[0], flux_neighbour_coord_tupel[1]])
-                    # due to the warning in the documentation of np.allclose
-                    # that np.allclose is not symetric, we test if both
-                    # np.allclose(a,b) and np.allclose(b,a) are true.
-                    a_close_b = np.allclose(flux_dof_coord, flux_neighbour_coord, rtol=1e-10, atol=1e-14)
-                    b_close_a = np.allclose(flux_neighbour_coord, flux_dof_coord, rtol=1e-10, atol=1e-14)
-                    if a_close_b and b_close_a:
-                        gravity_neumann_flux_neighbour_dof = flux_neighbour_dof
-                        break
-
-                if gravity_neumann_flux_neighbour_dof is None:
-                    raise RuntimeError(f"didn't find gravity_neumann_flux_neighbour_dof corresponding to dict key ({flux_dof_coord})",
-                                        "something is wrong")
-
-                gravity_neumann_flux.vector()[flux_dof_index] = gravity_neumann_flux_neighbour_dof
-
-            # In the following case, the gravity_neumann_flux of the neighbour is saved to
-            # gravity_neumann_flux currently and will be overwritten in the next step,
-            # if we don't save it to gravity_neumann_flux_prev of the neighbour.
-            if iteration == neighbour_iter_num:
-                # print(f"we are on subdomain{subdomain} and interface{interface.global_index} (alt index:{ind}) has neighbouring \n",
-                #       f"subdomains: subdomain{subdomain} and subdomain{neighbour}")
-                interface.gravity_neumann_flux_prev[neighbour][phase].update(#
-                    {global_facet_index: interface.gravity_neumann_flux[neighbour][phase][global_facet_index].copy()}
-                    )
-
-        return gravity_neumann_flux
-    ### END read_gravity_neumann_flux
-
-
     def calc_gl0_term(self, debug: bool = False) -> None: # tp.Dict[str, fl.Form]
         """calculate the gl0 terms for all interfaces of the subdomain.
 
@@ -671,10 +457,10 @@ class DomainPatch(df.SubDomain):
             S = self.saturation
 
             phases_for_gl0 = self.has_phases
-            # if self.isRichards and self.TPR_occures:
-            #     interface_has_TP_neighbour = not interface.neighbour_isRichards[self.subdomain_index]
-            #     if interface_has_TP_neighbour and self.include_gravity:
-            #         phases_for_gl0 = ['wetting', 'nonwetting']
+            if self.isRichards and self.TPR_occures:
+                interface_has_TP_neighbour = not interface.neighbour_isRichards[self.subdomain_index]
+                if interface_has_TP_neighbour and self.include_gravity:
+                    phases_for_gl0 = ['wetting', 'nonwetting']
 
             for phase in phases_for_gl0:
 
@@ -683,17 +469,18 @@ class DomainPatch(df.SubDomain):
                 # relative permeability
                 ka = self.relative_permeability[phase]
                 # previous timestep pressure
-                # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
-                #     pa_prev = 0
-                # else:
-                pa_prev = self.pressure_prev_timestep[phase]
+                if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+                    # dummy definition. we don't use this in this case
+                    pa_prev = 0
+                else:
+                    pa_prev = self.pressure_prev_timestep[phase]
 
                 if self.include_gravity:
                     z_a = self.gravity_term[phase]
+
                 if phase == 'wetting':
                     # the wetting phase is always present and we always need
                     # to calculate a gl0 term.
-                    # TODO: add intrinsic permeabilty!!!!
                     if self.include_gravity:
                         # z_a included the minus sign, see self._calc_gravity_expressions
                         flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev - z_a)
@@ -701,17 +488,16 @@ class DomainPatch(df.SubDomain):
                         flux = -1/mu_a*ka(S(pc_prev))*df.grad(pa_prev)
                 else:
                     # phase == 'nonwetting'
-                    # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
-                    #     if self.include_gravity:
-                    #         flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
-                    #     else:
-                    #         flux = 0
-                    # else:
-                    # we need anothre flux in this case
-                    if self.include_gravity:
-                        flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
+                    if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting":
+                        if self.include_gravity:
+                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(-z_a)
+                        else:
+                            flux = 0
                     else:
-                        flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
+                        if self.include_gravity:
+                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev - z_a)
+                        else:
+                            flux = -1/mu_a*ka(1-S(pc_prev))*df.grad(pa_prev)
 
                 flux_function = df.project(flux, self.function_space["flux"][phase])
                 flux_x, flux_y = flux_function.split()
@@ -738,15 +524,15 @@ class DomainPatch(df.SubDomain):
                     for gli_dof_index, gli_dof_coord in gli_dofs_along_facet.items():
                         flux_x_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_x"]
                         flux_y_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["flux_y"]
-                        p_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["pressure"]
-
-                        # if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
-                        #     gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
-                        #         + flux_y.vector()[flux_y_dof_index]*normal[1]
-                        # else:
-                        gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
-                            + flux_y.vector()[flux_y_dof_index]*normal[1] \
-                            - Lambda[phase]*pa_prev.vector()[p_dof_index]
+
+                        if self.isRichards and interface_has_TP_neighbour and phase is "nonwetting" and self.include_gravity:
+                            gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+                                + flux_y.vector()[flux_y_dof_index]*normal[1]
+                        else:
+                            p_dof_index = corresponding_dof_index[local_facet_ind][gli_dof_index]["pressure"]
+                            gl0.vector()[gli_dof_index] = flux_x.vector()[flux_x_dof_index]*normal[0] \
+                                + flux_y.vector()[flux_y_dof_index]*normal[1] \
+                                - Lambda[phase]*pa_prev.vector()[p_dof_index]
                         if debug:
                             print(f"gl0.vector()[gli_dof_ind] = {gl0.vector()[gli_dof_index]}")
 
@@ -1356,12 +1142,12 @@ class DomainPatch(df.SubDomain):
             marker_value = self.interface[interf].marker_value
             space = self.function_space
 
-            neighbour_assumes_R = self.interface[interf].neighbour_isRichards[self.subdomain_index]
+            neighbour_assumes_TP = not self.interface[interf].neighbour_isRichards[self.subdomain_index]
             # if our domain assumes Richards but the neighbour of the interface
             # assumes TP, we need to initialise interface values for the nonwetting
             # phase, too for communication, see the article.
-            if self.isRichards and not neighbour_assumes_R:
-                subdom_has_phases=["wetting","nonwetting"]
+            if self.isRichards and neighbour_assumes_TP:
+                subdom_has_phases=["wetting", "nonwetting"]
             else:
                 subdom_has_phases=self.has_phases
 
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
old mode 100644
new mode 100755
index 3fe92a0..13d24f6
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd-realistic.py
@@ -30,27 +30,26 @@ resolutions = {
                 # 4: 7e-7,
                 # 8: 7e-7,
                 # 16: 7e-7,
-                32: 1e-6,
+                32: 5e-6,
                 # 64: 7e-7,
                 # 128: 7e-7,
                 # 256: 7e-7
                 }
 
 ############ GRID #######################
-# mesh_resolution = 20
-timestep_size = 0.0005
-number_of_timesteps = 20
+timestep_size = 0.001
+number_of_timesteps = 15
 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 = 1
+number_of_timesteps_to_analyse = 5
 starttimes = [0.0]
 
-Lw = 10 #/timestep_size
-Lnw= 50
+Lw = 0.04 #/timestep_size
+Lnw= 0.8
 
-lambda_w = 40000
-lambda_nw = 40000
+lambda_w = 2
+lambda_nw = 8
 
 include_gravity = True
 debugflag = True
@@ -190,14 +189,21 @@ lambda_param = {#
          'nonwetting': lambda_nw}
 }
 
+intrinsic_permeability = {
+    1: {"wetting": 10e-6,
+        "nonwetting": 10e-6},
+    2: {"wetting": 10e-6,
+        "nonwetting": 10e-6},
+}
+
 ## relative permeabilty functions on subdomain 1
 def rel_perm1w(s):
     # relative permeabilty wetting on subdomain1
-    return s**2
+    return intrinsic_permeability[1]["wetting"]*s**2
 
 def rel_perm1nw(s):
     # relative permeabilty nonwetting on subdomain1
-    return (1-s)**2
+    return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
 
 _rel_perm1w = ft.partial(rel_perm1w)
 _rel_perm1nw = ft.partial(rel_perm1nw)
@@ -232,11 +238,11 @@ relative_permeability = {#
 # relative permeabilty functions on subdomain 1
 def rel_perm1w_prime(s):
     # relative permeabilty on subdomain1
-    return 2*s
+    return intrinsic_permeability[1]["wetting"]*2*s
 
 def rel_perm1nw_prime(s):
     # relative permeabilty on subdomain1
-    return -2*(1-s)
+    return -1*intrinsic_permeability[1]["nonwetting"]*2*(1-s)
 
 # # definition of the derivatives of the relative permeabilities
 # # relative permeabilty functions on subdomain 1
@@ -340,22 +346,22 @@ def saturation_sym_prime(pc, index):
 # note that the conditional definition of S-pc in the nonsymbolic part will be
 # incorporated in the construction of the exact solution below.
 S_pc_sym = {
-    1: ft.partial(saturation_sym, index=2),
-    2: ft.partial(saturation_sym, index=2),
+    1: ft.partial(saturation_sym, index=1),
+    2: ft.partial(saturation_sym, index=1),
     # 3: ft.partial(saturation_sym, index=2),
     # 4: ft.partial(saturation_sym, index=1)
 }
 
 S_pc_sym_prime = {
-    1: ft.partial(saturation_sym_prime, index=2),
-    2: ft.partial(saturation_sym_prime, index=2),
+    1: ft.partial(saturation_sym_prime, index=1),
+    2: ft.partial(saturation_sym_prime, index=1),
     # 3: ft.partial(saturation_sym_prime, index=2),
     # 4: ft.partial(saturation_sym_prime, index=1)
 }
 
 sat_pressure_relationship = {
-    1: ft.partial(saturation, index=2),
-    2: ft.partial(saturation, index=2),
+    1: ft.partial(saturation, index=1),
+    2: ft.partial(saturation, index=1),
     # 3: ft.partial(saturation, index=2),
     # 4: ft.partial(saturation, index=1)
 }
@@ -373,6 +379,7 @@ t = sym.symbols('t', positive=True)
 #         'nonwetting': (-1-t*(1.1+y + x**2))*y**3}, #*(1-x)**2*(1+x)**2*(1+y)**2},
 # } #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
 
+
 p_e_sym = {
     1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)),  #*cutoff,
         'nonwetting': (-1 -t*(1.1+ y*y))},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
diff --git a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
index bd55ce3..2c861d2 100755
--- a/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
+++ b/Two-phase-Richards/two-patch/TP-R-two-patch-test-case/TP-R-2-patch-pure-dd.py
@@ -30,30 +30,30 @@ resolutions = {
                 # 4: 7e-7,
                 # 8: 7e-7,
                 # 16: 7e-7,
-                32: 1e-6,
-                # 64: 7e-7,
+                # 32: 1e-6,
+                64: 7e-7,
                 # 128: 7e-7,
                 # 256: 7e-7
                 }
 
 ############ GRID #######################
 # mesh_resolution = 20
-timestep_size = 0.01
-number_of_timesteps = 150
+timestep_size = 0.005
+number_of_timesteps = 250
 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 = 5
 starttimes = [0.0]
 
-Lw = 0.1 #/timestep_size
-Lnw= 0.1
+Lw = 0.025 #/timestep_size
+Lnw= 0.025
 
-lambda_w = 4
-lambda_nw = 4
+lambda_w = 2
+lambda_nw = 2
 include_gravity = True
 debugflag = False
-analyse_condition = False
+analyse_condition = True
 
 output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
 
@@ -170,7 +170,7 @@ densities = {
         'nonwetting': 1} #1.225},
 }
 
-gravity_acceleration = 9.81
+gravity_acceleration = 1 #9.81
 
 L = {#
 # subdom_num : subdomain L for L-scheme
@@ -373,10 +373,10 @@ t = sym.symbols('t', positive=True)
 # } #-y*y*(sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)) - t*t*x*(0.5-y)*y*(1-x)
 
 p_e_sym = {
-    1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)),  #*cutoff,
-        'nonwetting': (-1 -t*(1.1+ y*y))},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
-    2: {'wetting': (-6.0 - (1.0 + t*t)*(1.0 + x*x)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
-        'nonwetting': (-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+    1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y))},  #*cutoff,
+        # 'nonwetting': (-1 -t*(1.1+ y*y))},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+    2: {'wetting': (-6 - (1+t*t)*(1 + x*x + y*y)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
+        'nonwetting': -t*(1.1 + y*x)*(sym.sin((x*y-0.5*t)*y**3))**4},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
 }
 
 pc_e_sym = dict()
diff --git a/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py b/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py
index 6206cd9..ef3b3d5 100755
--- a/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py
+++ b/Two-phase-Two-phase/two-patch/injection/TP-TP-2-patch-injection.py
@@ -38,7 +38,7 @@ resolutions = {
 
 ############ GRID #######################
 # mesh_resolution = 20
-timestep_size = 0.000001
+timestep_size = 0.0005
 number_of_timesteps = 13
 plot_timestep_every = 1
 # decide how many timesteps you want analysed. Analysed means, that we write out
@@ -46,11 +46,11 @@ plot_timestep_every = 1
 number_of_timesteps_to_analyse = 0
 starttime = 0.0
 
-Lw = 0.5 #/timestep_size
+Lw = 0.05 #/timestep_size
 Lnw=Lw
 
-lambda_w = 4000
-lambda_nw = 4000
+lambda_w = 80
+lambda_nw = 80
 
 include_gravity = False
 debugflag = True
@@ -177,6 +177,14 @@ densities = {
         'nonwetting': 1225}, #1225},
 }
 
+intrinsic_permeability = {
+    1: {"wetting": 10e-6,
+        "nonwetting": 10e-6},
+    2: {"wetting": 10e-6,
+        "nonwetting": 10e-6},
+}
+
+
 gravity_acceleration = 9.81
 
 L = {#
@@ -199,11 +207,11 @@ lambda_param = {#
 ## relative permeabilty functions on subdomain 1
 def rel_perm1w(s):
     # relative permeabilty wetting on subdomain1
-    return s**2
+    return intrinsic_permeability[1]["wetting"]*s**2
 
 def rel_perm1nw(s):
     # relative permeabilty nonwetting on subdomain1
-    return (1-s)**2
+    return intrinsic_permeability[1]["nonwetting"]*(1-s)**2
 
 _rel_perm1w = ft.partial(rel_perm1w)
 _rel_perm1nw = ft.partial(rel_perm1nw)
@@ -215,10 +223,10 @@ subdomain1_rel_perm = {
 ## relative permeabilty functions on subdomain 2
 def rel_perm2w(s):
     # relative permeabilty wetting on subdomain2
-    return s**3
+    return intrinsic_permeability[2]["wetting"]*s**3
 def rel_perm2nw(s):
     # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
-    return (1-s)**3
+    return intrinsic_permeability[2]["nonwetting"]*(1-s)**3
 
 _rel_perm2w = ft.partial(rel_perm2w)
 _rel_perm2nw = ft.partial(rel_perm2nw)
@@ -239,21 +247,21 @@ relative_permeability = {#
 # relative permeabilty functions on subdomain 1
 def rel_perm1w_prime(s):
     # relative permeabilty on subdomain1
-    return 2*s
+    return intrinsic_permeability[1]["wetting"]*2*s
 
 def rel_perm1nw_prime(s):
     # relative permeabilty on subdomain1
-    return -2*(1-s)
+    return -1*intrinsic_permeability[1]["nonwetting"]*2*(1-s)
 
 # # definition of the derivatives of the relative permeabilities
 # # relative permeabilty functions on subdomain 1
 def rel_perm2w_prime(s):
     # relative permeabilty on subdomain1
-    return 3*s**2
+    return intrinsic_permeability[2]["wetting"]*3*s**2
 
 def rel_perm2nw_prime(s):
     # relative permeabilty on subdomain1
-    return -3*(1-s)**2
+    return -3*intrinsic_permeability[2]["nonwetting"]*(1-s)**2
 
 _rel_perm1w_prime = ft.partial(rel_perm1w_prime)
 _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
@@ -378,8 +386,8 @@ t = sym.symbols('t', positive=True)
 initial_condition = {
     1: {'wetting': sym.printing.ccode(-6 - (1+t*t)*(1 + x*x + y*y)),  #*cutoff,
         'nonwetting': sym.printing.ccode(-1 -t*(1.1+ y*y))},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
-    2: {'wetting': sym.printing.ccode(-6.0 - (1.0 + t*t)*(1.0 + x*x)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
-        'nonwetting': sym.printing.ccode(-1 -t*(1.1 + y*y) - sym.sin((x*y-0.5*t)*y**2)**2)},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
+    2: {'wetting': sym.printing.ccode(-6 - (1+t*t)*(1 + x*x + y*y)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
+        'nonwetting': sym.printing.ccode(-1 -t*(1.1 + y*y))},  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2},
 }
 
 ### constructing source experessions.
@@ -408,7 +416,7 @@ extraction_radius = 0.075
 #
 def mollifier2d(x, y, epsilon):
     """ one d mollifier """
-    out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2))
+    out_expr = 0.05*sym.exp(-1/(1-(x**2 + y**2)/epsilon**2))
     return out_expr
 
 mollifier2d_handle_i = ft.partial(mollifier2d, epsilon=injection_radius)
-- 
GitLab