Skip to content
Snippets Groups Projects
Select Git revision
  • 5ca414d7b66c68cb46076b195432043a04533a3f
  • 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

RR-2-patch-test.py

Blame
  • RR-2-patch-test.py 13.99 KiB
    #!/usr/bin/python3
    import dolfin as df
    import mshr
    import numpy as np
    import sympy as sym
    import typing as tp
    import domainPatch as dp
    import LDDsimulation as ldd
    import functools as ft
    #import ufl as ufl
    
    ##### Domain and Interface ####
    # global simulation domain domain
    sub_domain0_vertices = [df.Point(0.0,0.0), #
                            df.Point(1.0,0.0),#
                            df.Point(1.0,1.0),#
                            df.Point(0.0,1.0)]
    # interface between subdomain1 and subdomain2
    interface12_vertices = [df.Point(0.0, 0.5),
                            df.Point(1.0, 0.5) ]
    # subdomain1.
    sub_domain1_vertices = [interface12_vertices[0],
                            interface12_vertices[1],
                            df.Point(1.0,1.0),
                            df.Point(0.0,1.0) ]
    # subdomain2
    sub_domain2_vertices = [df.Point(0.0,0.0),
                            df.Point(1.0,0.0),
                            interface12_vertices[1],
                            interface12_vertices[0] ]
    # list of subdomains given by the boundary polygon vertices.
    # Subdomains are given as a list of dolfin points forming
    # a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
    # to create the subdomain. subdomain_def_points[0] contains the
    # vertices of the global simulation domain and subdomain_def_points[i] contains the
    # vertices of the subdomain i.
    subdomain_def_points = [sub_domain0_vertices,#
                          sub_domain1_vertices,#
                          sub_domain2_vertices]
    # in the below list, index 0 corresponds to the 12 interface which has index 1
    interface_def_points = [interface12_vertices]
    # adjacent_subdomains[i] contains the indices of the subdomains sharing the
    # interface i (i.e. given by interface_def_points[i]).
    adjacent_subdomains = [[1,2]]
    isRichards = [True, True]
    
    Tmax = 2
    dt1 = 0.1
    dt2 = dt1
    # the timestep is taken as a dictionary to be able to run different timesteps on
    # different subdomains.
    timestep = {
        1: dt1,
        2: dt2
    }
    
    viscosity = {#
    # subdom_num : viscosity
        1 : {'wetting' :1}, #
        2 : {'wetting' :1}
    }
    
    porosity = {#
    # subdom_num : porosity
        1 : 1,#
        2 : 1
    }
    
    L = {#
    # subdom_num : subdomain L for L-scheme
        1 : {'wetting' :0.25},#
        2 : {'wetting' :0.25}
    }
    
    lambda_param = {#
    # subdom_num : lambda parameter for the L-scheme
        1 : {'wetting' :4},#
        2 : {'wetting' :4}
    }
    
    ## relative permeabilty functions on subdomain 1
    def rel_perm1(s):
        # relative permeabilty on subdomain1
        return s**2
    
    _rel_perm1 = ft.partial(rel_perm1)
    
    subdomain1_rel_perm = {
        'wetting': _rel_perm1,#
        'nonwetting': None
    }
    ## relative permeabilty functions on subdomain 2
    def rel_perm2(s):
        # relative permeabilty on subdomain2
        return s**3
    _rel_perm2 = ft.partial(rel_perm2)
    
    subdomain2_rel_perm = {
        'wetting': _rel_perm2,#
        'nonwetting': None
    }
    
    ## dictionary of relative permeabilties on all domains.
    relative_permeability = {#
        1: subdomain1_rel_perm,
        2: subdomain2_rel_perm
    }
    
    def saturation(pressure, subdomain_index):
        # inverse capillary pressure-saturation-relationship
        return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
    
    sat_pressure_relationship = {#
        1: ft.partial(saturation, subdomain_index = 1),#
        2: ft.partial(saturation, subdomain_index = 2)
    }
    
    
    
    # def saturation(pressure, subdomain_index):
    #     # inverse capillary pressure-saturation-relationship
    #     return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
    #
    # sa
    
    
    # initialise LDD simulation class
    simulation = ldd.LDDsimulation()
    simulation.set_parameters(output_dir = "",#
        subdomain_def_points = subdomain_def_points,#
        isRichards = isRichards,#
        interface_def_points = interface_def_points,#
        adjacent_subdomains = adjacent_subdomains,#
        mesh_resolution = 2,#
        viscosity = viscosity,#
        porosity = porosity,#
        L = L,#
        lambda_param = lambda_param,#
        relative_permeability = relative_permeability,#
        saturation = sat_pressure_relationship,#
        timestep = timestep#
        )
    
    simulation.initialise()
    # simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
    # subdomain marker functions
    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.set_all(0)
    # interface = dp.Interface(vertices=interface_def_points[0], #
    #                                     tol=mesh.hmin()/100,#
    #                                     internal=True,#
    #                                     adjacent_subdomains = adjacent_subdomains[0])
    # print("interface vertices are: \n", interface.vertices)
    # interface.mark(interface_marker, 1)
    # simulation._init_interfaces(interface_def_points, adjacent_subdomains)
    
    interface = simulation.interface
    interface_marker = simulation.interface_marker
    
    subdomain_boundary_marker1 = df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1)
    subdomain_boundary_marker2 = df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
    subdomain_boundary_marker1.set_all(0)
    subdomain_boundary_marker2.set_all(0)
    # print(dir(subdomain_boundary_marker1))
    
    interface[0].mark(subdomain_boundary_marker1, 1)
    interface[0].mark(subdomain_boundary_marker2, 1)
    
    
    # domain measures
    dx1 = df.Measure('dx', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1)
    dx2 = df.Measure('dx', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
    
    # boundary measures
    ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1)
    ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
    
    V1 = df.FunctionSpace(mesh_subdomain[1], 'P', 1)
    p1_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2, t = 0)
    
    coordinates = mesh_subdomain[0].coordinates()
    print("global mesh coordinates of domain: \n", coordinates)
    coordinates1 = mesh_subdomain[1].coordinates()
    
    coordinates2 = mesh_subdomain[2].coordinates()
    submesh1_data = mesh_subdomain[1].data()
    dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0)
    print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices)
    #
    # coordinates2 = mesh_subdomain[2].coordinates()
    # print("\nvertices of subdomain2: \n", coordinates2)
    submesh2_data = mesh_subdomain[2].data()
    dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
    # print("'parent_vertex_indices',1 \n")
    #
    #
    # print("on domain: \n")
    # interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
    # interface_coordinates = interface[0]._vertex_indices(interface_marker, 1, print_vertex_indices = True)
    print("\non subdomain2: \n")
    # dom2_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker2, 1, print_coordinates = True)
    dom2_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
    
    
    dofs_on_interface1 = interface[0].dofs_on_interface(V1, subdomain_boundary_marker1, 1)
    print("Dofs on Interface dom1", dofs_on_interface1)
    
    V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
    dofs_on_interface2 = interface[0].dofs_on_interface(V2, subdomain_boundary_marker2, 1)
    print("Dofs on Interface dom1", dofs_on_interface2)
    
    # markermesh = subdomain_boundary_marker2.mesh()
    # functionSpaceMesh = V2.mesh()
    # print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates())
    # print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates())
    
    dom1_d2v = df.dof_to_vertex_map(V1)
    dom1_v2d = df.vertex_to_dof_map(V1)
    print("dof to vertex map for V1:\n", dom1_d2v)
    print("vertex to dof map for V1:\n", dom1_v2d)
    dom2_d2v = df.dof_to_vertex_map(V2)
    dom2_v2d = df.vertex_to_dof_map(V2)
    print("dof to vertex map for V2:\n", dom2_d2v)
    print("vertex to dof map for V2:\n", dom2_v2d)
    
    testf1 = df.Constant(41)
    f1test = df.interpolate(testf1, V1)
    testf2 = df.Constant(42)
    f2test = df.interpolate(testf2, V2)
    
    f1vec = f1test.compute_vertex_values()
    for i in range(len(f1vec)):
        f1vec[i] = i
    
    # f1test.vector()[:] = f1vec
    f1nodal = f1test.vector().get_local()
    # for i in range(len(f1nodal)):
    #     f1nodal[i] = i
    
    print("nodal values of f1nodal\n", f1nodal)
    
    #
    # print('dofs on interface \n', f1vec[dofs_on_interface])
    # f1test.vector()[:] = f1vec
    # print("f1test.vector()[:]\n", f1test.vector()[:])
    # print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
    
    print("constant function on dom1\n", f1vec)
    
    for i in dofs_on_interface1:
        print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
        # f1test.vector() is enumerated by dof not by vertex
        f1test.vector()[i] = i
    # this enumerates according to the dof numbering
    print("testfunction f1test with set interface values\n", f1test.vector()[:])
    print("testfunction f1test with get_local\n", f1test.vector().get_local())
    # f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
    # .vector()[:] and .vector().get_local() order values according to the vertex numbering
    print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
    print("testfunction f1test with compute_vertex_values mit d2vmap\n", f1test.compute_vertex_values()[dom1_d2v[:]])
    i=0
    for x in coordinates1:
         print(f"vertex {i}: f1test[{dom1_v2d[i]}] = {f1test.vector()[dom1_v2d[i]]}\tu({x}) = {f1test(x)}")
         i += 1
    print("\non subdomain1: \n")
    print("vertices of subdomain1: \n", coordinates1)
    # dom1_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker1, 1, print_coordinates = True)
    dom1_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker1, 1, print_vertex_indices = True)
    
    interface[0].write_dofs(fem_function = f1test, #
                            interface_dofs = dofs_on_interface1,#
                            dof_to_vert_map = dom1_d2v,#
                            local_to_parent_vertex_map = dom1_parent_vertex_indices)
    
    interface[0].read_dofs(fem_function = f2test, #
                            interface_dofs = dofs_on_interface2,#
                            dof_to_vert_map = dom2_d2v,#
                            local_to_parent_vertex_map = dom2_parent_vertex_indices)
    
    i=0
    for x in coordinates2:
         print(f"vertex {i}: f2test[{dom2_v2d[i]}] = {f2test.vector()[dom2_v2d[i]]}\tu({x}) = {f2test(x)}")
         i += 1
    print("\non subdomain2: \n")
    print("vertices of subdomain2: \n", coordinates2)
    
    # dom2_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker2, 1, print_coordinates = True)
    dom2_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
    print("Parent indices", dom2_parent_vertex_indices)
    
    
    u1 = df.TrialFunction(V1)
    v1 = df.TestFunction(V1)
    
    u2 = df.TrialFunction(V2)
    v2 = df.TestFunction(V2)
    
    def relative_permeability(s, subdomain_index):
        if subdomain_index == 1:
            # relative permeabilty on subdomain1
            return s**2
        else:
            # relative permeabilty on subdomain2
            return s**3
    
    def saturation(pressure, subdomain_index):
        # inverse capillary pressure-saturation-relationship
        return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1)
    
    
    # exact solution
    # cell = df.Cell("triangle", 2)
    # x1 = df.SpatialCoordinate(mesh_subdomain[1])
    # x2 = df.SpatialCoordinate(mesh_subdomain[2])
    
    p2_exact = df.Expression('1.0 - (1.0 + t*t)*(1.0 + x[1]*x[1])', domain = mesh_subdomain[2], degree = 2, t = 0)
    p1_e = df.interpolate(p1_exact, V1)
    p2_e = df.interpolate(p2_exact, V2)
    #
    # initial conditions
    p1_initial = df.Expression('-(x[0]*x[0] + x[1]*x[1])', domain = mesh_subdomain[1], degree = 2)
    p2_initial = df.Expression('-x[1]*x[1]', domain = mesh_subdomain[2], degree = 2)
    p1_0 = df.interpolate(p1_initial, V1)
    p2_0 = df.interpolate(p2_initial, V2)
    
    # Initial boundary value for the pressure on interface
    p_gamma1 = df.DirichletBC(V1, df.Constant(0.5), subdomain_boundary_marker1, 1)
    p_gamma2 = df.DirichletBC(V2, df.Constant(0.1), subdomain_boundary_marker2, 1)
    
    outerBC1 = df.DirichletBC(V1, p1_exact, subdomain_boundary_marker1, 0)
    outerBC2 = df.DirichletBC(V2, p2_exact, subdomain_boundary_marker2, 0)
    
    
    ### source terms
    source1 = df.Expression('4.0/pow(1 + x[0]*x[0] + x[1]*x[1], 2) - t/sqrt( pow(1 + t*t, 3)*(1 + x[0]*x[0] + x[1]*x[1]) )', domain = mesh_subdomain[1], degree = 2, t = 0)
    source2 = df.Expression('2.0*(1-x[1]*x[1])/pow(1 + x[1]*x[1], 2) - 2*t/(3*pow( pow(1 + t*t, 4)*(1 + x[1]*x[1]), 1/3))', domain = mesh_subdomain[2], degree = 2, t = 0)
    source1h = df.interpolate(source1, V1)
    source2h = df.interpolate(source2, V2)
    
    L1 = 0.25
    L2 = L1
    timestep = 0.1
    lambda1 = 4
    lambda2 = lambda1
    # initial iteartion
    p1_i = p1_0
    p2_i = p2_0
    # initial g_i
    g1_i = -lambda1*p1_i
    g2_i = -lambda1*p2_i
    
    for iterations in range(1):
        a1 = (L1*u1*v1)*dx1 + (timestep*df.dot(relative_permeability(saturation(p1_i, 1), 1)*df.grad(u1), df.grad(v1)))*dx1 + (timestep*lambda1*u1*v1)*ds1
        rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 +(timestep*source1*v1)*dx1 -(timestep*(g1_i)*v1)*ds1
        A1 = df.assemble(a1)
        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('solution:\n', U1)
        p1_i.vector()[:] = U1
        interface[0].write_dofs(fem_function = p1_i, #
                                interface_dofs = dofs_on_interface1,#
                                dof_to_vert_map = dom1_d2v,#
                                local_to_parent_vertex_map = dom1_parent_vertex_indices)
    #
    # Save mesh to file
    df.File('./test_domain_layered_soil.xml.gz') << 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_boundary_marker.pvd') << subdomain_boundary_marker1
    df.File('./test_subdomain2_boundary_marker.pvd') << subdomain_boundary_marker2