Skip to content
Snippets Groups Projects
Commit 5b53e577 authored by David Seus's avatar David Seus
Browse files

just the current iteration of the RR-test-file

parent 3a7fc0d5
Branches
Tags
No related merge requests found
...@@ -6,6 +6,7 @@ import sympy as sym ...@@ -6,6 +6,7 @@ import sympy as sym
import typing as tp import typing as tp
import domainPatch as dp import domainPatch as dp
import LDDsimulation as ldd import LDDsimulation as ldd
import functools as ft
#import ufl as ufl #import ufl as ufl
##### Domain and Interface #### ##### Domain and Interface ####
...@@ -29,22 +30,83 @@ sub_domain2_vertices = [df.Point(0.0,0.0), ...@@ -29,22 +30,83 @@ sub_domain2_vertices = [df.Point(0.0,0.0),
interface12_vertices[0] ] interface12_vertices[0] ]
# list of subdomains given by the boundary polygon vertices. # list of subdomains given by the boundary polygon vertices.
# Subdomains are given as a list of dolfin points forming # Subdomains are given as a list of dolfin points forming
# a closed polygon, such that mshr.Polygon(subdomain_vertices[i]) can be used # a closed polygon, such that mshr.Polygon(subdomain_def_points[i]) can be used
# to create the subdomain. subdomain_vertices[0] contains the # to create the subdomain. subdomain_def_points[0] contains the
# vertices of the global simulation domain and subdomain_vertices[i] contains the # vertices of the global simulation domain and subdomain_def_points[i] contains the
# vertices of the subdomain i. # vertices of the subdomain i.
subdomain_vertices = [sub_domain0_vertices,# subdomain_def_points = [sub_domain0_vertices,#
sub_domain1_vertices,# sub_domain1_vertices,#
sub_domain2_vertices] sub_domain2_vertices]
# in the below list, index 0 corresponds to the 12 interface which has index 1 # in the below list, index 0 corresponds to the 12 interface which has index 1
interface_vertices = [interface12_vertices] interface_def_points = [interface12_vertices]
# adjacent_subdomains[i] contains the indices of the subdomains sharing the # adjacent_subdomains[i] contains the indices of the subdomains sharing the
# interface i (i.e. given by interface_vertices[i]). # interface i (i.e. given by interface_def_points[i]).
adjacent_subdomains = [[1,2]] adjacent_subdomains = [[1,2]]
is_Richards = [True, True]
viscosity = {#
# subdom_num : viscosity
1 : [1], #
2 : [1]
}
porosity = {#
# subdom_num : porosity
1 : 1,#
2 : 1
}
L = {#
# subdom_num : subdomain L for L-scheme
1 : [0.25],#
2 : [0.25]
}
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
1 : [4],#
2 : [4]
}
def rel_perm1(s):
# relative permeabilty on subdomain1
return s**2
_rel_perm1 = ft.partial(rel_perm1)
def rel_perm2(s):
# relative permeabilty on subdomain2
return s**3
_rel_perm2 = ft.partial(rel_perm2)
relative_permeability = {#
1: _rel_perm1,
2: _rel_perm2,
}
# 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 # initialise LDD simulation class
simulation = ldd.LDDsimulation() simulation = ldd.LDDsimulation()
simulation.init_meshes_and_markers(subdomain_vertices, mesh_resolution=2) simulation.set_parameters(output_dir = "",#
subdomain_def_points = subdomain_def_points,#
is_Richards = is_Richards,#
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)
simulation.initialise()
# simulation._init_meshes_and_markers(subdomain_def_points, mesh_resolution=2)
# subdomain marker functions # subdomain marker functions
domain_marker = simulation.domain_marker domain_marker = simulation.domain_marker
mesh_subdomain = simulation.mesh_subdomain mesh_subdomain = simulation.mesh_subdomain
...@@ -52,232 +114,232 @@ mesh_subdomain = simulation.mesh_subdomain ...@@ -52,232 +114,232 @@ mesh_subdomain = simulation.mesh_subdomain
# mesh = mesh_subdomain[0] # mesh = mesh_subdomain[0]
# interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1) # interface_marker = df.MeshFunction('size_t', mesh, mesh.topology().dim()-1)
# interface_marker.set_all(0) # interface_marker.set_all(0)
# interface = dp.Interface(vertices=interface_vertices[0], # # interface = dp.Interface(vertices=interface_def_points[0], #
# tol=mesh.hmin()/100,# # tol=mesh.hmin()/100,#
# internal=True,# # internal=True,#
# adjacent_subdomains = adjacent_subdomains[0]) # adjacent_subdomains = adjacent_subdomains[0])
# print("interface vertices are: \n", interface.vertices) # print("interface vertices are: \n", interface.vertices)
# interface.mark(interface_marker, 1) # interface.mark(interface_marker, 1)
simulation.init_interfaces(interface_vertices, adjacent_subdomains) # simulation._init_interfaces(interface_def_points, adjacent_subdomains)
interface = simulation.interface interface = simulation.interface
interface_marker = simulation.interface_marker interface_marker = simulation.interface_marker
#
# subdomain_boundary_marker1 = df.MeshFunction('size_t', mesh_subdomain[1], mesh_subdomain[1].topology().dim()-1) 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_marker2 = df.MeshFunction('size_t', mesh_subdomain[2], mesh_subdomain[2].topology().dim()-1)
# subdomain_boundary_marker1.set_all(0) subdomain_boundary_marker1.set_all(0)
# subdomain_boundary_marker2.set_all(0) subdomain_boundary_marker2.set_all(0)
# # print(dir(subdomain_boundary_marker1)) # print(dir(subdomain_boundary_marker1))
#
# interface[0].mark(subdomain_boundary_marker1, 1) interface[0].mark(subdomain_boundary_marker1, 1)
# interface[0].mark(subdomain_boundary_marker2, 1) interface[0].mark(subdomain_boundary_marker2, 1)
#
#
# # domain measures # domain measures
# dx1 = df.Measure('dx', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1) 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) dx2 = df.Measure('dx', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
#
# # boundary measures # boundary measures
# ds1 = df.Measure('ds', domain = mesh_subdomain[1], subdomain_data = subdomain_boundary_marker1) 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) ds2 = df.Measure('ds', domain = mesh_subdomain[2], subdomain_data = subdomain_boundary_marker2)
#
# V1 = df.FunctionSpace(mesh_subdomain[1], 'P', 1) 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) 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() coordinates = mesh_subdomain[0].coordinates()
# print("global mesh coordinates of domain: \n", coordinates) print("global mesh coordinates of domain: \n", coordinates)
# coordinates1 = mesh_subdomain[1].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() # coordinates2 = mesh_subdomain[2].coordinates()
# submesh1_data = mesh_subdomain[1].data() # print("\nvertices of subdomain2: \n", coordinates2)
# dom1_parent_vertex_indices = submesh1_data.array('parent_vertex_indices',0) submesh2_data = mesh_subdomain[2].data()
# print("map of dom1 vertices to parent vertices on global domain:\n", dom1_parent_vertex_indices) dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0)
# # # print("'parent_vertex_indices',1 \n")
# # coordinates2 = mesh_subdomain[2].coordinates() #
# # print("\nvertices of subdomain2: \n", coordinates2) #
# submesh2_data = mesh_subdomain[2].data() # print("on domain: \n")
# dom2_parent_vertex_indices = submesh2_data.array('parent_vertex_indices',0) # interface_coordinates = interface[0].coordinates(interface_marker, 1, print_coordinates = True)
# # print("'parent_vertex_indices',1 \n") # 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)
# # print("on domain: \n") dom2_interface_def_points = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
# # 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") dofs_on_interface1 = interface[0].dofs_on_interface(V1, subdomain_boundary_marker1, 1)
# # dom2_interface_coordinates = interface[0].coordinates(subdomain_boundary_marker2, 1, print_coordinates = True) print("Dofs on Interface dom1", dofs_on_interface1)
# dom2_interface_vertices = interface[0]._vertex_indices(subdomain_boundary_marker2, 1, print_vertex_indices = True)
# V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1)
# dofs_on_interface2 = interface[0].dofs_on_interface(V2, subdomain_boundary_marker2, 1)
# dofs_on_interface1 = interface[0].dofs_on_interface(V1, subdomain_boundary_marker1, 1) print("Dofs on Interface dom1", dofs_on_interface2)
# print("Dofs on Interface dom1", dofs_on_interface1)
# # markermesh = subdomain_boundary_marker2.mesh()
# V2 = df.FunctionSpace(mesh_subdomain[2], 'P', 1) # functionSpaceMesh = V2.mesh()
# dofs_on_interface2 = interface[0].dofs_on_interface(V2, subdomain_boundary_marker2, 1) # print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates())
# print("Dofs on Interface dom1", dofs_on_interface2) # print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates())
#
# # markermesh = subdomain_boundary_marker2.mesh() dom1_d2v = df.dof_to_vertex_map(V1)
# # functionSpaceMesh = V2.mesh() dom1_v2d = df.vertex_to_dof_map(V1)
# # print("Domain 2 mesh extracted from markerfunction:\n", markermesh.coordinates()) print("dof to vertex map for V1:\n", dom1_d2v)
# # print("Domain 2 mesh extracted from Function Space V2:\n", functionSpaceMesh.coordinates()) print("vertex to dof map for V1:\n", dom1_v2d)
# dom2_d2v = df.dof_to_vertex_map(V2)
# dom1_d2v = df.dof_to_vertex_map(V1) dom2_v2d = df.vertex_to_dof_map(V2)
# dom1_v2d = df.vertex_to_dof_map(V1) print("dof to vertex map for V2:\n", dom2_d2v)
# print("dof to vertex map for V1:\n", dom1_d2v) print("vertex to dof map for V2:\n", dom2_v2d)
# print("vertex to dof map for V1:\n", dom1_v2d)
# dom2_d2v = df.dof_to_vertex_map(V2) testf1 = df.Constant(41)
# dom2_v2d = df.vertex_to_dof_map(V2) f1test = df.interpolate(testf1, V1)
# print("dof to vertex map for V2:\n", dom2_d2v) testf2 = df.Constant(42)
# print("vertex to dof map for V2:\n", dom2_v2d) f2test = df.interpolate(testf2, V2)
#
# testf1 = df.Constant(41) f1vec = f1test.compute_vertex_values()
# f1test = df.interpolate(testf1, V1) for i in range(len(f1vec)):
# testf2 = df.Constant(42) f1vec[i] = i
# f2test = df.interpolate(testf2, V2)
# # f1test.vector()[:] = f1vec
# f1vec = f1test.compute_vertex_values() f1nodal = f1test.vector().get_local()
# for i in range(len(f1vec)): # for i in range(len(f1nodal)):
# f1vec[i] = i # f1nodal[i] = i
#
# # f1test.vector()[:] = f1vec print("nodal values of f1nodal\n", f1nodal)
# 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_vertices = 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_vertices = 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): # print('dofs on interface \n', f1vec[dofs_on_interface])
# 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 # f1test.vector()[:] = f1vec
# rhs1 = (L1*p1_i*v1)*dx1 - ((saturation(p1_i, 1) - saturation(p1_0, 1))*v1)*dx1 + (timestep*(source1 - g1_i)*v1)*ds1 # print("f1test.vector()[:]\n", f1test.vector()[:])
# A1 = df.assemble(a1) # print("f1test.vector()[dofs_on_interface]\n", f1test.vector()[dofs_on_interface])
# b1 = df.assemble(rhs1)
# p_gamma1.apply(A1,b1) print("constant function on dom1\n", f1vec)
# outerBC1.apply(A1,b1)
# u = df.Function(V1) for i in dofs_on_interface1:
# U1 = u.vector() print("local dof ind %i, parent = %i", i,dom1_parent_vertex_indices[dom1_d2v[i]])
# df.solve(A1,U1,b1) # f1test.vector() is enumerated by dof not by vertex
# print('solution:\n', U1) f1test.vector()[i] = i
# p1_i.vector()[:] = U1 # this enumerates according to the dof numbering
# interface[0].write_dofs(fem_function = p1_i, # print("testfunction f1test with set interface values\n", f1test.vector()[:])
# interface_dofs = dofs_on_interface1,# print("testfunction f1test with get_local\n", f1test.vector().get_local())
# dof_to_vert_map = dom1_d2v,# # f1test.vector()[dom1_v2d[dom1_d2v[dofs_on_interface]]] = f1vec[dofs_on_interface]
# local_to_parent_vertex_map = dom1_parent_vertex_indices) # .vector()[:] and .vector().get_local() order values according to the vertex numbering
# # print("testfunction f1test with compute_vertex_values\n", f1test.compute_vertex_values())
# # Save mesh to file 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 - 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_domain_layered_soil.xml.gz') << mesh_subdomain[0]
df.File('./test_global_interface_marker.pvd') << interface_marker df.File('./test_global_interface_marker.pvd') << interface_marker
# df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1] df.File('./test_subdomain1.xml.gz') << mesh_subdomain[1]
# df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2] df.File('./test_subdomain2.xml.gz') << mesh_subdomain[2]
# df.File('./test_domain_marker.pvd') << domain_marker df.File('./test_domain_marker.pvd') << domain_marker
# df.File('./test_domain_layered_soil_solution.pvd') << u df.File('./test_domain_layered_soil_solution.pvd') << u
# df.File('./test_subdomain1_boundary_marker.pvd') << subdomain_boundary_marker1 df.File('./test_subdomain1_boundary_marker.pvd') << subdomain_boundary_marker1
# df.File('./test_subdomain2_boundary_marker.pvd') << subdomain_boundary_marker2 df.File('./test_subdomain2_boundary_marker.pvd') << subdomain_boundary_marker2
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment