Select Git revision
RR-2-patch-test.py
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