Skip to content
Snippets Groups Projects
Select Git revision
  • dde2e70da0d92e96a7e3cbd789dcda29e419ecad
  • master default protected
  • parallelistation_of_subdomain_loop
  • parallelistation_test
  • nonzero_gli_term_for_nonwetting_in_TPR
  • testing_updating_pwi_as_pwiminus1_into_calculation_of_pnwi
  • v2.2
  • v2.1
  • v2.0
  • v1.0
10 results

TP-TP-layered_soil_with_inner_patch.py

Blame
  • TP-TP-layered_soil_with_inner_patch.py 25.08 KiB
    #!/usr/bin/python3
    """This program sets up a domain together with a decomposition into subdomains
    modelling layered soil. This is used for our LDD article with tp-tp and tp-r
    coupling.
    
    Along with the subdomains and the mesh domain markers are set upself.
    The resulting mesh is saved into files for later use.
    """
    
    #!/usr/bin/python3
    import dolfin as df
    import mshr
    import numpy as np
    import sympy as sym
    import typing as tp
    import functools as ft
    import domainPatch as dp
    import LDDsimulation as ldd
    import helpers as hlp
    
    # init sympy session
    sym.init_printing()
    
    use_case = "TP-TP-layered-soil-with-inner-patch"
    solver_tol = 1E-6
    
    ############ GRID #######################ü
    mesh_resolution = 50
    timestep_size = 0.001
    number_of_timesteps = 1500
    # 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 = 10
    starttime = 0
    
    Lw = 0.25  #/timestep_size
    Lnw=Lw
    
    lambda_w = 41
    lambda_nw = 41
    
    include_gravity = True
    debugflag = False
    analyse_condition = False
    
    output_string = "./output/nondirichlet_number_of_timesteps{}_".format(number_of_timesteps)
    
    # global domain
    subdomain0_vertices = [df.Point(-1.0,-1.0), #
                            df.Point(1.0,-1.0),#
                            df.Point(1.0,1.0),#
                            df.Point(-1.0,1.0)]
    
    interface12_vertices = [df.Point(-1.0, 0.8),
                            df.Point(0.3, 0.8),
                            df.Point(0.5, 0.9),
                            df.Point(0.8, 0.7),
                            df.Point(1.0, 0.65)]
    
    
                            # interface23
    interface23_vertices = [df.Point(-1.0, 0.0),
                            df.Point(-0.35, 0.0),
                            # df.Point(6.5, 4.5),
                            df.Point(0.0, 0.0)]
    
    interface24_vertices = [interface23_vertices[2],
                            df.Point(0.6, 0.0),
                            ]
    
    interface25_vertices = [interface24_vertices[1],
                            df.Point(1.0, 0.0)
                            ]
    
    
    interface32_vertices = [interface23_vertices[2],
                            interface23_vertices[1],
                            interface23_vertices[0]]
    
    
    interface36_vertices = [df.Point(-1.0, -0.6),
                            df.Point(-0.6, -0.45)]
    
    
    interface46_vertices = [interface36_vertices[1],
                            df.Point(0.3, -0.25)]
    
    interface56_vertices = [interface46_vertices[1],
                            df.Point(0.65, -0.6),
                            df.Point(1.0, -0.7)]
    
    
    
    
    interface34_vertices = [interface36_vertices[1],
                            interface23_vertices[2]]
    # interface36
    
    interface45_vertices = [interface56_vertices[0],
                            df.Point(0.7, -0.2),
                            interface25_vertices[0]
                            ]
    
    # # subdomain1.
    # subdomain1_vertices = [interface12_vertices[0],
    #                         interface12_vertices[1],
    #                         interface12_vertices[2],
    #                         interface12_vertices[3],
    #                         interface12_vertices[4], # southern boundary, 12 interface
    #                         subdomain0_vertices[2], # eastern boundary, outer boundary
    #                         subdomain0_vertices[3]] # northern boundary, outer on_boundary
    #
    # # vertex coordinates of the outer boundaries. If it can not be specified as a
    # # polygon, use an entry per boundary polygon. This information is used for defining
    # # the Dirichlet boundary conditions. If a domain is completely internal, the
    # # dictionary entry should be 0: None
    # subdomain1_outer_boundary_verts = {
    #     0: [interface12_vertices[4], #
    #         subdomain0_vertices[2], # eastern boundary, outer boundary
    #         subdomain0_vertices[3],
    #         interface12_vertices[0]]
    # }
    #
    
    
    # #subdomain1
    # subdomain2_vertices = [interface23_vertices[0],
    #                         interface23_vertices[1],
    #                         interface23_vertices[2],
    #                         interface23_vertices[3],
    #                         interface23_vertices[4],
    #                         interface23_vertices[5], # southern boundary, 23 interface
    #                         subdomain1_vertices[4], # eastern boundary, outer boundary
    #                         subdomain1_vertices[3],
    #                         subdomain1_vertices[2],
    #                         subdomain1_vertices[1],
    #                         subdomain1_vertices[0] ] # northern boundary, 12 interface
    #
    # subdomain2_outer_boundary_verts = {
    #     0: [interface23_vertices[5],
    #         subdomain1_vertices[4]],
    #     1: [subdomain1_vertices[0],
    #         interface23_vertices[0]]
    # }
    #
    
    # interface_vertices introduces a global numbering of interfaces.
    interface_def_points = [interface12_vertices,
                            interface23_vertices,
                            interface24_vertices,
                            interface25_vertices,
                            interface34_vertices,
                            interface36_vertices,
                            interface45_vertices,
                            interface46_vertices,
                            interface56_vertices,
                            ]
    adjacent_subdomains = [[1,2],
                           [2,3],
                           [2,4],
                           [2,5],
                           [3,4],
                           [3,6],
                           [4,5],
                           [4,6],
                           [5,6]
                           ]
    
    # subdomain1.
    subdomain1_vertices = [interface12_vertices[0],
                            interface12_vertices[1],
                            interface12_vertices[2],
                            interface12_vertices[3],
                            interface12_vertices[4], # southern boundary, 12 interface
                            subdomain0_vertices[2], # eastern boundary, outer boundary
                            subdomain0_vertices[3]] # northern boundary, outer on_boundary
    
    # vertex coordinates of the outer boundaries. If it can not be specified as a
    # polygon, use an entry per boundary polygon. This information is used for defining
    # the Dirichlet boundary conditions. If a domain is completely internal, the
    # dictionary entry should be 0: None
    subdomain1_outer_boundary_verts = {
        0: [subdomain1_vertices[4], #
            subdomain1_vertices[5], # eastern boundary, outer boundary
            subdomain1_vertices[6],
            subdomain1_vertices[0]]
    }
    
    #subdomain1
    subdomain2_vertices = [interface23_vertices[0],
                            interface23_vertices[1],
                            interface23_vertices[2],
                            interface24_vertices[1],
                            interface25_vertices[1], # southern boundary, 23 interface
                            subdomain1_vertices[4], # eastern boundary, outer boundary
                            subdomain1_vertices[3],
                            subdomain1_vertices[2],
                            subdomain1_vertices[1],
                            subdomain1_vertices[0] ] # northern boundary, 12 interface
    
    subdomain2_outer_boundary_verts = {
        0: [subdomain2_vertices[9],
            subdomain2_vertices[0]],
        1: [subdomain2_vertices[4],
            subdomain2_vertices[5]]
    }
    
    
    subdomain3_vertices = [interface36_vertices[0],
                           interface36_vertices[1],
                           # interface34_vertices[0],
                           interface34_vertices[1],
                           # interface32_vertices[0],
                           interface32_vertices[1],
                           interface32_vertices[2]
                           ]
    
    subdomain3_outer_boundary_verts = {
        0: [subdomain3_vertices[4],
            subdomain3_vertices[0]]
    }
    
    
    # subdomain3
    subdomain4_vertices = [interface46_vertices[0],
                           interface46_vertices[1],
                           interface45_vertices[1],
                           interface24_vertices[1],
                           interface24_vertices[0],
                           interface34_vertices[1]
                           ]
    
    subdomain4_outer_boundary_verts = None
    
    subdomain5_vertices = [interface56_vertices[0],
                           interface56_vertices[1],
                           interface56_vertices[2],
                           interface25_vertices[1],
                           interface25_vertices[0],
                           interface45_vertices[1],
                           interface45_vertices[0]
    ]
    
    subdomain5_outer_boundary_verts = {
        0: [subdomain5_vertices[2],
            subdomain5_vertices[3]]
    }
    
    
    
    subdomain6_vertices = [subdomain0_vertices[0],
                           subdomain0_vertices[1], # southern boundary, outer boundary
                           interface56_vertices[2],
                           interface56_vertices[1],
                           interface56_vertices[0],
                           interface36_vertices[1],
                           interface36_vertices[0]
                           ]
    
    subdomain6_outer_boundary_verts = {
        0: [subdomain6_vertices[6],
            subdomain6_vertices[0],
            subdomain6_vertices[1],
            subdomain6_vertices[2]]
    }
    
    
    subdomain_def_points = [subdomain0_vertices,#
                          subdomain1_vertices,#
                          subdomain2_vertices,#
                          subdomain3_vertices,#
                          subdomain4_vertices,
                          subdomain5_vertices,
                          subdomain6_vertices
                          ]
    
    
    # if a subdomain has no outer boundary write None instead, i.e.
    # i: None
    # if i is the index of the inner subdomain.
    outer_boundary_def_points = {
        # subdomain number
        1: subdomain1_outer_boundary_verts,
        2: subdomain2_outer_boundary_verts,
        3: subdomain3_outer_boundary_verts,
        4: subdomain4_outer_boundary_verts,
        5: subdomain5_outer_boundary_verts,
        6: subdomain6_outer_boundary_verts
    }
    
    
    isRichards = {
        1: False,
        2: False,
        3: False,
        4: False,
        5: False,
        6: False
        }
    
    # isRichards = {
    #     1: True,
    #     2: True,
    #     3: True,
    #     4: True,
    #     5: True,
    #     6: True
    #     }
    
    # Dict of the form: { subdom_num : viscosity }
    viscosity = {
        1: {'wetting' :1,
             'nonwetting': 1/50},
        2: {'wetting' :1,
             'nonwetting': 1/50},
        3: {'wetting' :1,
             'nonwetting': 1/50},
        4: {'wetting' :1,
             'nonwetting': 1/50},
        5: {'wetting' :1,
             'nonwetting': 1/50},
        6: {'wetting' :1,
             'nonwetting': 1/50},
    }
    
    # Dict of the form: { subdom_num : density }
    densities = {
        1: {'wetting': 1,  #997
             'nonwetting': 1},  #1},  #1.225},
        2: {'wetting': 1,  #997
             'nonwetting': 1},  #1.225},
        3: {'wetting': 1,  #997
             'nonwetting': 1},  #1.225},
        4: {'wetting': 1,  #997
             'nonwetting': 1},  #1.225}
        5: {'wetting': 1,  #997
             'nonwetting': 1},  #1.225},
        6: {'wetting': 1,  #997
             'nonwetting': 1}  #1.225}
    }
    
    gravity_acceleration = 9.81
    # porosities taken from
    # https://www.geotechdata.info/parameter/soil-porosity.html
    # Dict of the form: { subdom_num : porosity }
    porosity = {
        1: 1,  #0.2,  # Clayey gravels, clayey sandy gravels
        2: 1,  #0.22, # Silty gravels, silty sandy gravels
        3: 1,  #0.37, # Clayey sands
        4: 1,  #0.2 # Silty or sandy clay
        5: 1,  #
        6: 1,  #
    }
    
    # subdom_num : subdomain L for L-scheme
    L = {
        1: {'wetting' :Lw,
             'nonwetting': Lnw},
        2: {'wetting' :Lw,
             'nonwetting': Lnw},
        3: {'wetting' :Lw,
             'nonwetting': Lnw},
        4: {'wetting' :Lw,
             'nonwetting': Lnw},
        5: {'wetting' :Lw,
             'nonwetting': Lnw},
        6: {'wetting' :Lw,
             'nonwetting': Lnw}
    }
    
    # subdom_num : lambda parameter for the L-scheme
    lambda_param = {
        1: {'wetting': lambda_w,
             'nonwetting': lambda_nw},#
        2: {'wetting': lambda_w,
             'nonwetting': lambda_nw},#
        3: {'wetting': lambda_w,
             'nonwetting': lambda_nw},#
        4: {'wetting': lambda_w,
             'nonwetting': lambda_nw},#
        5: {'wetting': lambda_w,
             'nonwetting': lambda_nw},#
        6: {'wetting': lambda_w,
             'nonwetting': lambda_nw},#
    }
    
    
    ## relative permeabilty functions on subdomain 1
    def rel_perm1w(s):
        # relative permeabilty wetting on subdomain1
        return s**2
    
    
    def rel_perm1nw(s):
        # relative permeabilty nonwetting on subdomain1
        return (1-s)**2
    
    
    ## relative permeabilty functions on subdomain 2
    def rel_perm2w(s):
        # relative permeabilty wetting on subdomain2
        return 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
    
    
    _rel_perm1w = ft.partial(rel_perm1w)
    _rel_perm1nw = ft.partial(rel_perm1nw)
    _rel_perm2w = ft.partial(rel_perm2w)
    _rel_perm2nw = ft.partial(rel_perm2nw)
    
    subdomain1_rel_perm = {
        'wetting': _rel_perm1w,#
        'nonwetting': _rel_perm1nw
    }
    
    subdomain2_rel_perm = {
        'wetting': _rel_perm2w,#
        'nonwetting': _rel_perm2nw
    }
    
    # _rel_perm3 = ft.partial(rel_perm2)
    # subdomain3_rel_perm = subdomain2_rel_perm.copy()
    #
    # _rel_perm4 = ft.partial(rel_perm1)
    # subdomain4_rel_perm = subdomain1_rel_perm.copy()
    
    # dictionary of relative permeabilties on all domains.
    relative_permeability = {
        1: subdomain1_rel_perm,
        2: subdomain1_rel_perm,
        3: subdomain2_rel_perm,
        4: subdomain2_rel_perm,
        5: subdomain2_rel_perm,
        6: subdomain2_rel_perm,
    }
    
    # definition of the derivatives of the relative permeabilities
    # relative permeabilty functions on subdomain 1
    def rel_perm1w_prime(s):
        # relative permeabilty on subdomain1
        return 2*s
    
    def rel_perm1nw_prime(s):
        # relative permeabilty on subdomain1
        return -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
    
    def rel_perm2nw_prime(s):
        # relative permeabilty on subdomain1
        return -3*(1-s)**2
    
    _rel_perm1w_prime = ft.partial(rel_perm1w_prime)
    _rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
    _rel_perm2w_prime = ft.partial(rel_perm2w_prime)
    _rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)
    
    subdomain1_rel_perm_prime = {
        'wetting': _rel_perm1w_prime,
        'nonwetting': _rel_perm1nw_prime
    }
    
    
    subdomain2_rel_perm_prime = {
        'wetting': _rel_perm2w_prime,
        'nonwetting': _rel_perm2nw_prime
    }
    
    # dictionary of relative permeabilties on all domains.
    ka_prime = {
        1: subdomain1_rel_perm_prime,
        2: subdomain1_rel_perm_prime,
        3: subdomain2_rel_perm_prime,
        4: subdomain2_rel_perm_prime,
        5: subdomain2_rel_perm_prime,
        6: subdomain2_rel_perm_prime,
    }
    
    
    
    # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
    # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
    # this function needs to be monotonically decreasing in the capillary pressure pc.
    # since in the richards case pc=-pw, this becomes as a function of pw a mono
    # tonically INCREASING function like in our Richards-Richards paper. However
    # since we unify the treatment in the code for Richards and two-phase, we need
    # the same requierment
    # for both cases, two-phase and Richards.
    # def saturation(pc, n_index, alpha):
    #     # inverse capillary pressure-saturation-relationship
    #     return df.conditional(pc > 0, 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index)), 1)
    #
    # # S-pc-relation ship. We use the van Genuchten approach, i.e. pc = 1/alpha*(S^{-1/m} -1)^1/n, where
    # # we set alpha = 0, assume m = 1-1/n (see Helmig) and assume that residual saturation is Sw
    # def saturation_sym(pc, n_index, alpha):
    #     # inverse capillary pressure-saturation-relationship
    #     #df.conditional(pc > 0,
    #     return 1/((1 + (alpha*pc)**n_index)**((n_index - 1)/n_index))
    #
    #
    # # derivative of S-pc relationship with respect to pc. This is needed for the
    # # construction of a analytic solution.
    # def saturation_sym_prime(pc, n_index, alpha):
    #     # inverse capillary pressure-saturation-relationship
    #     return -(alpha*(n_index - 1)*(alpha*pc)**(n_index - 1)) / ( (1 + (alpha*pc)**n_index)**((2*n_index - 1)/n_index) )
    #
    # derivative of S-pc relationship with respect to pc. This is needed for the
    # construction of a analytic solution.
    
    #
    # # 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, n_index=3, alpha=0.001),
    #     2: ft.partial(saturation_sym, n_index=3, alpha=0.001),
    #     3: ft.partial(saturation_sym, n_index=3, alpha=0.001),
    #     4: ft.partial(saturation_sym, n_index=3, alpha=0.001),
    #     5: ft.partial(saturation_sym, n_index=3, alpha=0.001),
    #     6: ft.partial(saturation_sym, n_index=3, alpha=0.001)
    # }
    #
    # S_pc_sym_prime = {
    #     1: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
    #     2: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
    #     3: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
    #     4: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
    #     5: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001),
    #     6: ft.partial(saturation_sym_prime, n_index=3, alpha=0.001)
    # }
    #
    # sat_pressure_relationship = {
    #     1: ft.partial(saturation, n_index=3, alpha=0.001),
    #     2: ft.partial(saturation, n_index=3, alpha=0.001),
    #     3: ft.partial(saturation, n_index=3, alpha=0.001),
    #     4: ft.partial(saturation, n_index=3, alpha=0.001),
    #     5: ft.partial(saturation, n_index=3, alpha=0.001),
    #     6: ft.partial(saturation, n_index=3, alpha=0.001)
    # }
    
    def saturation(pc, n_index):
        # inverse capillary pressure-saturation-relationship
        return df.conditional(pc > 0, 1/((1 + pc)**(1/(n_index + 1))), 1)
    
    
    def saturation_sym(pc, n_index):
        # inverse capillary pressure-saturation-relationship
        return 1/((1 + pc)**(1/(n_index + 1)))
    
    def saturation_sym_prime(pc, n_index):
        # inverse capillary pressure-saturation-relationship
        return -1/((n_index+1)*(1 + pc)**((n_index+2)/(n_index+1)))
    
    
    S_pc_sym = {
        1: ft.partial(saturation_sym, n_index=1),
        2: ft.partial(saturation_sym, n_index=1),
        3: ft.partial(saturation_sym, n_index=2),
        4: ft.partial(saturation_sym, n_index=2),
        5: ft.partial(saturation_sym, n_index=2),
        6: ft.partial(saturation_sym, n_index=2)
    }
    
    S_pc_sym_prime = {
        1: ft.partial(saturation_sym_prime, n_index=1),
        2: ft.partial(saturation_sym_prime, n_index=1),
        3: ft.partial(saturation_sym_prime, n_index=2),
        4: ft.partial(saturation_sym_prime, n_index=2),
        5: ft.partial(saturation_sym_prime, n_index=2),
        6: ft.partial(saturation_sym_prime, n_index=2)
    }
    
    sat_pressure_relationship = {
        1: ft.partial(saturation, n_index=1),
        2: ft.partial(saturation, n_index=1),
        3: ft.partial(saturation, n_index=2),
        4: ft.partial(saturation, n_index=2),
        5: ft.partial(saturation, n_index=2),
        6: ft.partial(saturation, n_index=2)
    }
    
    
    #############################################
    # Manufacture source expressions with sympy #
    #############################################
    x, y = sym.symbols('x[0], x[1]')  # needed by UFL
    t = sym.symbols('t', positive=True)
    
    
    p_e_sym = {
        1: {'wetting': -5.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
            'nonwetting': (-1 -t*(1.1 + y + x**2)) },
        2: {'wetting': -5.0 - (1.0 + t*t)*(1.0 + x*x + y*y),
            'nonwetting': (-1 -t*(1.1 + y + x**2)) },
        3: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
            'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
        4: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
            'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
        5: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
            'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
        6: {'wetting': (-5.0 - (1.0 + t*t)*(1.0 + x*x)),
            'nonwetting': (-1 -t*(1 + x**2) - sym.sqrt(2+t**2)*(1+y)*y**2) },
        # 2: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)),
        #     'nonwetting': - 2 - t*(1 + (y-5.0) + x**2)**2 -sym.sqrt(2+t**2)*(1 + (y-5.0))},
        # 3: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)*3*sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)),
        #     'nonwetting': - 2 - t*(1 + x**2)**2 -sym.sqrt(2+t**2)},
        # 4: {'wetting': 1.0 - (1.0 + t*t)*(10.0 + x*x + (y-5.0)*(y-5.0)*3*sym.sin(-2*t+2*x)*sym.sin(1/2*y-1.2*t)),
        #     'nonwetting': - 2 - t*(1 + x**2)**2 -sym.sqrt(2+t**2)}
    }
    
    # pc_e_sym = {
    #     1: p_e_sym[1]['nonwetting'] - p_e_sym[1]['wetting'],
    #     2: p_e_sym[2]['nonwetting'] - p_e_sym[2]['wetting'],
    #     3: p_e_sym[3]['nonwetting'] - p_e_sym[3]['wetting'],
    #     4: p_e_sym[4]['nonwetting'] - p_e_sym[4]['wetting'],
    #     5: p_e_sym[5]['nonwetting'] - p_e_sym[5]['wetting'],
    #     6: p_e_sym[5]['nonwetting'] - p_e_sym[6]['wetting']
    # }
    
    # pc_e_sym = {
    #     1: -p_e_sym[1]['wetting'],
    #     2: -p_e_sym[2]['wetting'],
    #     3: -p_e_sym[3]['wetting'],
    #     4: -p_e_sym[4]['wetting'],
    #     5: -p_e_sym[5]['wetting'],
    #     6: -p_e_sym[6]['wetting']
    # }
    
    
    pc_e_sym = dict()
    for subdomain, isR in isRichards.items():
        if isR:
            pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
        else:
            pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
                                            - p_e_sym[subdomain]['wetting']})
    
    
    symbols = {"x": x,
               "y": y,
               "t": t}
    # turn above symbolic code into exact solution for dolphin and
    # construct the rhs that matches the above exact solution.
    exact_solution_example = hlp.generate_exact_solution_expressions(
                            symbols=symbols,
                            isRichards=isRichards,
                            symbolic_pressure=p_e_sym,
                            symbolic_capillary_pressure=pc_e_sym,
                            saturation_pressure_relationship=S_pc_sym,
                            saturation_pressure_relationship_prime=S_pc_sym_prime,
                            viscosity=viscosity,
                            porosity=porosity,
                            relative_permeability=relative_permeability,
                            relative_permeability_prime=ka_prime,
                            densities=densities,
                            gravity_acceleration=gravity_acceleration,
                            include_gravity=include_gravity,
                            )
    source_expression = exact_solution_example['source']
    exact_solution = exact_solution_example['exact_solution']
    initial_condition = exact_solution_example['initial_condition']
    
    # Dictionary of dirichlet boundary conditions.
    dirichletBC = dict()
    # similarly to the outer boundary dictionary, if a patch has no outer boundary
    # None should be written instead of an expression.
    # This is a bit of a brainfuck:
    # dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
    # Since a domain patch can have several disjoint outer boundary parts, the
    # expressions need to get an enumaration index which starts at 0.
    # So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
    # subdomain ind and boundary part j.
    # Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
    # return the actual expression needed for the dirichlet condition for both
    # phases if present.
    
    # subdomain index: {outer boudary part index: {phase: expression}}
    for subdomain in isRichards.keys():
        # if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None
        if outer_boundary_def_points[subdomain] is None:
            dirichletBC.update({subdomain: None})
        else:
            dirichletBC.update({subdomain: dict()})
            # set the dirichlet conditions to be the same code as exact solution on
            # the subdomain.
            for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
                dirichletBC[subdomain].update(
                    {outer_boundary_ind: exact_solution[subdomain]}
                    )
    
    write_to_file = {
        'meshes_and_markers': True,
        'L_iterations': True
    }
    
    # initialise LDD simulation class
    simulation = ldd.LDDsimulation(tol=1E-14, debug=debugflag, LDDsolver_tol=solver_tol)
    simulation.set_parameters(use_case = use_case,
                              output_dir=output_string,
                              subdomain_def_points=subdomain_def_points,
                              isRichards=isRichards,
                              interface_def_points=interface_def_points,
                              outer_boundary_def_points=outer_boundary_def_points,
                              adjacent_subdomains=adjacent_subdomains,
                              mesh_resolution=mesh_resolution,
                              viscosity=viscosity,
                              porosity=porosity,
                              L=L,
                              lambda_param=lambda_param,
                              relative_permeability=relative_permeability,
                              saturation=sat_pressure_relationship,
                              starttime=starttime,
                              number_of_timesteps=number_of_timesteps,
                              number_of_timesteps_to_analyse=number_of_timesteps_to_analyse,
                              timestep_size=timestep_size,
                              sources=source_expression,
                              initial_conditions=initial_condition,
                              dirichletBC_expression_strings=dirichletBC,
                              exact_solution=exact_solution,
                              densities=densities,
                              include_gravity=include_gravity,
                              write2file=write_to_file,
                              )
    
    simulation.initialise()
    # print(simulation.__dict__)
    simulation.run(analyse_condition=analyse_condition)
    # simulation.LDDsolver(time=0, debug=True, analyse_timestep=True)
    # df.info(parameters, True)