diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py index 96681d1d5e12c039c913ad40e97ed2914541b1ba..2d676444e3508eda92f2f06c6e96e45b1db4c188 100644 --- a/LDDsimulation/LDDsimulation.py +++ b/LDDsimulation/LDDsimulation.py @@ -302,7 +302,7 @@ class LDDsimulation(object): # solution case. self._init_exact_solution_expression() self._init_DirichletBC_dictionary() - if self.exact_solution: + if self.exact_solution and self.write2file['solutions']: df.info(colored("writing exact solution for all time steps to xdmf files ...", "yellow")) self.write_exact_solution_to_xdmf() self._initialised = True @@ -379,7 +379,7 @@ class LDDsimulation(object): # for phase in subdomain.has_phases: # print(f"calculated pressure of {phase}phase:") # print(subdomain.pressure[phase].vector()[:]) - if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True): + if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True) and self.write2file['solutions']: subdomain.write_solution_to_xdmf(# file = self.solution_file[subdom_ind], # time = self.t,# @@ -417,7 +417,7 @@ class LDDsimulation(object): subdomain_index = subdom_ind, errors = relative_L2_errornorm, ) - if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True): + if np.isin(self.timestep_num, self.timesteps_to_plot, assume_unique=True) and self.write2file['absolute_differences']: pressure_exact.update( {phase: df.interpolate(pa_exact, subdomain.function_space["pressure"][phase])} ) @@ -455,6 +455,7 @@ class LDDsimulation(object): # df.info(colored("start post processing calculations ...\n", "yellow")) # self.post_processing() # df.info(colored("finished post processing calculations \nAll right. I'm Done.", "green")) + return self.spacetime_errornorm def LDDsolver(self, time: float = None, # debug: bool = False, # @@ -529,16 +530,17 @@ class LDDsimulation(object): if analyse_timestep: # save the calculated L2 errors to csv file for plotting # with pgfplots - subsequent_error_filename = self.output_dir\ - +self.output_filename_parameter_part[sd_index]\ - +"subsequent_iteration_errors" +"_at_time"+\ - "{number}".format(number=self.timestep_num) +".csv" #"{number:.{digits}f}".format(number=time, digits=4) - self.write_subsequent_errors_to_csv( - filename = subsequent_error_filename, # - subdomain_index = sd_index, - errors = subsequent_iter_error - ) - if analyse_condition: + if self.write2file['subsequent_errors']: + subsequent_error_filename = self.output_dir\ + +self.output_filename_parameter_part[sd_index]\ + +"subsequent_iteration_errors" +"_at_time"+\ + "{number}".format(number=self.timestep_num) +".csv" #"{number:.{digits}f}".format(number=time, digits=4) + self.write_subsequent_errors_to_csv( + filename = subsequent_error_filename, # + subdomain_index = sd_index, + errors = subsequent_iter_error + ) + if analyse_condition and self.write2file['condition_numbers']: # save the calculated condition numbers of the assembled # matrices a separate file for monitoring condition_number_filename = self.output_dir\ diff --git a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-horizontal-interface-avoiding-origin.py b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-horizontal-interface-avoiding-origin.py new file mode 100755 index 0000000000000000000000000000000000000000..e3a5a8e9c079712c14aa4ea75d958350e91fbfdd --- /dev/null +++ b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-horizontal-interface-avoiding-origin.py @@ -0,0 +1,556 @@ +#!/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 helpers as hlp +import datetime + +date = datetime.datetime.now() +datestr = date.strftime("%Y-%m-%d") +#import ufl as ufl + +# init sympy session +sym.init_printing() + +use_case = "TP-TP-2-patch-really-pure-dd-horizontal-interface-avoding-origin" +solver_tol = 1E-6 +max_iter_num = 1000 +FEM_Lagrange_degree = 1 + +############ GRID ####################### +mesh_resolution = 10 +timestep_size = 0.0001 +number_of_timesteps = 4000 +# 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 = 6 +starttime = 0 + +Lw = 0.25 #/timestep_size +Lnw=Lw + +lambda_w = 4 +lambda_nw = 4 + +include_gravity = False +debugflag = False +analyse_condition = True + +output_string = "./output/{}-{}_timesteps{}_".format(datestr, use_case, number_of_timesteps) + +##### Domain and Interface #### +# global simulation domain domain +sub_domain0_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)] +# interface between subdomain1 and subdomain2 +interface12_vertices = [df.Point(-1.0, 0.6), + df.Point(1.0, 0.6) ] +# subdomain1. +sub_domain1_vertices = [interface12_vertices[0], + interface12_vertices[1], + sub_domain0_vertices[2], + sub_domain0_vertices[3] ] + +# 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[1], + sub_domain0_vertices[2], + sub_domain0_vertices[3], # + interface12_vertices[0]] +} +# subdomain2 +sub_domain2_vertices = [sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1], + interface12_vertices[0] ] + +subdomain2_outer_boundary_verts = { + 0: [interface12_vertices[0], # + sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1]] +} +# subdomain2_outer_boundary_verts = { +# 0: [interface12_vertices[0], df.Point(0.0,0.0)],# +# 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], # +# 2: [df.Point(1.0,0.0), interface12_vertices[1]] +# } +# subdomain2_outer_boundary_verts = { +# 0: None +# } + +# 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] + +# 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 +} + +# 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 = { + 1: False, # + 2: False + } + + +viscosity = {# +# subdom_num : viscosity + 1 : {'wetting' :1, + 'nonwetting': 1}, # + 2 : {'wetting' :1, + 'nonwetting': 1} +} + +porosity = {# +# subdom_num : porosity + 1 : 1,# + 2 : 1 +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 1, #997, + 'nonwetting': 1}, #1225}, + 2: {'wetting': 1, #997, + 'nonwetting': 1}, #1225}, +} + +gravity_acceleration = 9.81 + + +L = {# +# subdom_num : subdomain L for L-scheme + 1 : {'wetting' :Lw, + 'nonwetting': Lnw},# + 2 : {'wetting' :Lw, + 'nonwetting': Lnw} +} + + +lambda_param = {# +# subdom_num : lambda parameter for the L-scheme + 1 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw},# + 2 : {'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 + +_rel_perm1w = ft.partial(rel_perm1w) +_rel_perm1nw = ft.partial(rel_perm1nw) + +subdomain1_rel_perm = { + 'wetting': _rel_perm1w,# + 'nonwetting': _rel_perm1nw +} +## relative permeabilty functions on subdomain 2 +def rel_perm2w(s): + # relative permeabilty wetting on subdomain2 + return s**2 +def rel_perm2nw(s): + # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2 + return (1-s)**2 + +_rel_perm2w = ft.partial(rel_perm2w) +_rel_perm2nw = ft.partial(rel_perm2nw) + +subdomain2_rel_perm = { + 'wetting': _rel_perm2w,# + 'nonwetting': _rel_perm2nw +} + +## dictionary of relative permeabilties on all domains. +relative_permeability = {# + 1: subdomain1_rel_perm, + 2: 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 2*s + +def rel_perm2nw_prime(s): + # relative permeabilty on subdomain1 + return -2*(1-s) + +_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: subdomain2_rel_perm_prime, +} + + + +def saturation(pc, index): + # inverse capillary pressure-saturation-relationship + return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1) + + + +def pc_sat_rel_sym(S, index): + # capillary pressure-saturation-relationship + return 1/S**(index+1) -1 + +pc_saturation_sym = { + 1: ft.partial(pc_sat_rel_sym, index=1), + 2: ft.partial(pc_sat_rel_sym, index=1), +} + +def saturation_sym(pc, index): + # inverse capillary pressure-saturation-relationship + return 1/((1 + pc)**(1/(index + 1))) + + +# derivative of S-pc relationship with respect to pc. This is needed for the +# construction of a analytic solution. +def saturation_sym_prime(pc, index): + # inverse capillary pressure-saturation-relationship + return -1/((index+1)*(1 + pc)**((index+2)/(index+1))) + + +# 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=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=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=1), + 2: ft.partial(saturation, index=1), + # 3: ft.partial(saturation, index=2), + # 4: ft.partial(saturation, index=1) +} + +# +# 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) ) +# +# # 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=6, 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=6, 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=6, alpha=0.001),p1w + Spc[1] +# # 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) +# } +# + + +############################################# +# Manufacture source expressions with sympy # +############################################# +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +symbols = { "x": x, + "y": y, + "t": t} + +# epsilon_x_inner = 0.7 +# epsilon_x_outer = 0.99 +# epsilon_y_inner = epsilon_x_inner +# epsilon_y_outer = epsilon_x_outer +# +# def mollifier(x, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1) +# return out_expr +# +# mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner) +# +# pw_sym_x = sym.Piecewise( +# (mollifier_handle(x), x**2 < epsilon_x_outer**2), +# (0, True) +# ) +# pw_sym_y = sym.Piecewise( +# (mollifier_handle(y), y**2 < epsilon_y_outer**2), +# (0, True(-1 -t*(1.1 + x**2) - sym.sqrt(2+t**2)*(1.1+y)**2*(0.6-y)**2)) +# ) +# +# def mollifier2d(x, y, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1) +# return out_expr +# +# mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer) +# +# pw_sym2d_x = sym.Piecewise( +# (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2), +# (0, True) +# ) +# +# zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise( +# (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))), +# (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise( +# (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))), +# (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise( +# (1, y<=-2*epsilon_x_inner), +# (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))), +# (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y +# gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x +# cutoff = gaussian/(gaussian + zero_on_shrinking) +# +# +# sat_sym = { +# 1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t), +# 2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t) +# } +# +# Spc = { +# 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)), +# 2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True)) +# } +# +# p1w = (-1 - (1+t*t)*(1 + x*x + y*y))#*cutoff +# p2w = p1w +# p_e_sym = { +# 1: {'wetting': p1w, +# 'nonwetting': (p1w + Spc[1])}, #*cutoff}, +# 2: {'wetting': p2w, +# 'nonwetting': (p2w + Spc[2])}, #*cutoff}, +# } + +p_e_sym = { + 1: {'wetting': (-6 - (1+t*t)*(1 + x*x + (0.6-y)**2)), #*cutoff, + 'nonwetting': (-1 -t*(1.1+ 0.6-y + x**2))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, + 2: {'wetting': (-6 - (1+t*t)*(1 + x*x + (0.6-y)**2)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2, + 'nonwetting': (-1 -t*(1.1+ 0.6-y + x**2))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, + # 1: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + # 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff}, + # 2: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + # 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff}, +} + + +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']}) + + + +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]} + ) + + +# def saturation(pressure, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) +# +# sa + +write_to_file = { + 'meshes_and_markers': True, + 'L_iterations': True +} + + +# initialise LDD simulation class +simulation = ldd.LDDsimulation( + tol=1E-14, + LDDsolver_tol=solver_tol, + debug=debugflag, + max_iter_num=max_iter_num, + FEM_Lagrange_degree=FEM_Lagrange_degree + ) + +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() +# simulation.write_exact_solution_to_xdmf() +simulation.run(analyse_condition=analyse_condition) diff --git a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-tilted-interface/TP-TP-2-patch-pure-dd-horizontal-tilted-interface-avoiding-origin.py b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-tilted-interface/TP-TP-2-patch-pure-dd-horizontal-tilted-interface-avoiding-origin.py new file mode 100755 index 0000000000000000000000000000000000000000..89d01b9a57ea283544de539e2b2e53611f68852f --- /dev/null +++ b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/TP-TP-2-patch-pure-dd-tilted-interface/TP-TP-2-patch-pure-dd-horizontal-tilted-interface-avoiding-origin.py @@ -0,0 +1,555 @@ +#!/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 helpers as hlp +import datetime + +date = datetime.datetime.now() +datestr = date.strftime("%Y-%m-%d") +#import ufl as ufl + +# init sympy session +sym.init_printing() + +use_case = "Actual-pure-dd-Tilted-interface-avoding-origin-same-space-for-pwpnw-nonwetting-first" +solver_tol = 2E-6 +max_iter_num = 1000 + +############ GRID ####################### +mesh_resolution = 20 +timestep_size = 0.0005 +number_of_timesteps = 1000 +# 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 = 6 +starttime = 0 + +Lw = 2 #/timestep_size +Lnw=Lw + +lambda_w = 4 +lambda_nw = 4 + +include_gravity = False +debugflag = False +analyse_condition = True + +output_string = "./output/{}-{}_timesteps{}_".format(datestr, use_case, number_of_timesteps) + +##### Domain and Interface #### +# global simulation domain domain +sub_domain0_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)] +# interface between subdomain1 and subdomain2 +interface12_vertices = [df.Point(-1.0, 0.9), + df.Point(1.0, 0.4) ] + # interface equation: y = -1/4*x + 13/20 +# subdomain1. +sub_domain1_vertices = [interface12_vertices[0], + interface12_vertices[1], + sub_domain0_vertices[2], + sub_domain0_vertices[3] ] + +# 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[1], + sub_domain0_vertices[2], + sub_domain0_vertices[3], # + interface12_vertices[0]] +} +# subdomain2 +sub_domain2_vertices = [sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1], + interface12_vertices[0] ] + +subdomain2_outer_boundary_verts = { + 0: [interface12_vertices[0], # + sub_domain0_vertices[0], + sub_domain0_vertices[1], + interface12_vertices[1]] +} +# subdomain2_outer_boundary_verts = { +# 0: [interface12_vertices[0], df.Point(0.0,0.0)],# +# 1: [df.Point(0.0,0.0), df.Point(1.0,0.0)], # +# 2: [df.Point(1.0,0.0), interface12_vertices[1]] +# } +# subdomain2_outer_boundary_verts = { +# 0: None +# } + +# 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] + +# 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 +} + +# 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 = { + 1: False, # + 2: False + } + + +viscosity = {# +# subdom_num : viscosity + 1 : {'wetting' :1, + 'nonwetting': 1}, # + 2 : {'wetting' :1, + 'nonwetting': 1} +} + +porosity = {# +# subdom_num : porosity + 1 : 1,# + 2 : 1 +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting': 1, #997, + 'nonwetting': 1}, #1225}, + 2: {'wetting': 1, #997, + 'nonwetting': 1}, #1225}, +} + +gravity_acceleration = 9.81 + + +L = {# +# subdom_num : subdomain L for L-scheme + 1 : {'wetting' :Lw, + 'nonwetting': Lnw},# + 2 : {'wetting' :Lw, + 'nonwetting': Lnw} +} + + +lambda_param = {# +# subdom_num : lambda parameter for the L-scheme + 1 : {'wetting' :lambda_w, + 'nonwetting': lambda_nw},# + 2 : {'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 + +_rel_perm1w = ft.partial(rel_perm1w) +_rel_perm1nw = ft.partial(rel_perm1nw) + +subdomain1_rel_perm = { + 'wetting': _rel_perm1w,# + 'nonwetting': _rel_perm1nw +} +## relative permeabilty functions on subdomain 2 +def rel_perm2w(s): + # relative permeabilty wetting on subdomain2 + return s**2 +def rel_perm2nw(s): + # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2 + return (1-s)**2 + +_rel_perm2w = ft.partial(rel_perm2w) +_rel_perm2nw = ft.partial(rel_perm2nw) + +subdomain2_rel_perm = { + 'wetting': _rel_perm2w,# + 'nonwetting': _rel_perm2nw +} + +## dictionary of relative permeabilties on all domains. +relative_permeability = {# + 1: subdomain1_rel_perm, + 2: 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 2*s + +def rel_perm2nw_prime(s): + # relative permeabilty on subdomain1 + return -2*(1-s) + +_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: subdomain2_rel_perm_prime, +} + + + +def saturation(pc, index): + # inverse capillary pressure-saturation-relationship + return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1) + + + +def pc_sat_rel_sym(S, index): + # capillary pressure-saturation-relationship + return 1/S**(index+1) -1 + +pc_saturation_sym = { + 1: ft.partial(pc_sat_rel_sym, index=1), + 2: ft.partial(pc_sat_rel_sym, index=1), +} + +def saturation_sym(pc, index): + # inverse capillary pressure-saturation-relationship + return 1/((1 + pc)**(1/(index + 1))) + + +# derivative of S-pc relationship with respect to pc. This is needed for the +# construction of a analytic solution. +def saturation_sym_prime(pc, index): + # inverse capillary pressure-saturation-relationship + return -1/((index+1)*(1 + pc)**((index+2)/(index+1))) + + +# 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=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=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=1), + 2: ft.partial(saturation, index=1), + # 3: ft.partial(saturation, index=2), + # 4: ft.partial(saturation, index=1) +} + +# +# 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) ) +# +# # 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=6, 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=6, 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=6, alpha=0.001),p1w + Spc[1] +# # 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) +# } +# + + +############################################# +# Manufacture source expressions with sympy # +############################################# +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +symbols = { "x": x, + "y": y, + "t": t} + +# epsilon_x_inner = 0.7 +# epsilon_x_outer = 0.99 +# epsilon_y_inner = epsilon_x_inner +# epsilon_y_outer = epsilon_x_outer +# +# def mollifier(x, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1) +# return out_expr +# +# mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner) +# +# pw_sym_x = sym.Piecewise( +# (mollifier_handle(x), x**2 < epsilon_x_outer**2), +# (0, True) +# ) +# pw_sym_y = sym.Piecewise( +# (mollifier_handle(y), y**2 < epsilon_y_outer**2), +# (0, True) +# ) +# +# def mollifier2d(x, y, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1) +# return out_expr +# +# mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer) +# +# pw_sym2d_x = sym.Piecewise( +# (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2), +# (0, True) +# ) +# +# zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise( +# (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))), +# (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise( +# (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))), +# (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise( +# (1, y<=-2*epsilon_x_inner), +# (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))), +# (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y +# gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x +# cutoff = gaussian/(gaussian + zero_on_shrinking) +# +# +# sat_sym = { +# 1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t), +# 2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t) +# } +# +# Spc = { +# 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)), +# 2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True)) +# } +# +# p1w = (-1 - (1+t*t)*(1 + x*x + y*y))#*cutoff +# p2w = p1w +# p_e_sym = { +# 1: {'wetting': p1w, +# 'nonwetting': (p1w + Spc[1])}, #*cutoff}, +# 2: {'wetting': p2w, +# 'nonwetting': (p2w + Spc[2])}, #*cutoff}, +# } + +p_e_sym = { + 1: {'wetting': (-6 - (1+t*t)*(1 + x*x + (0.6-y)**2)), #*cutoff, + 'nonwetting': (-1 -t*(1.1+ 0.6-y + x**2))}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, + 2: {'wetting': (-6 - (1+t*t)*(1 + x*x + (0.6-y)**2)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2, + 'nonwetting': (-1 -t*(1.1+ 0.6-y + x**2))}, + # 1: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + # 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff}, + # 2: {'wetting': (-5 - (1+t*t)*(1 + x*x + y*y)), #*cutoff, + # 'nonwetting': (-1 -t*(1.1+y + x**2))}, #*cutoff}, +} + + +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']}) + + + +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]} + ) + + +# def saturation(pressure, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) +# +# sa + +write_to_file = { + 'meshes_and_markers': True, + 'L_iterations': True +} + + +# initialise LDD simulation class +simulation = ldd.LDDsimulation( + tol=1E-14, + LDDsolver_tol=solver_tol, + debug=debugflag, + max_iter_num=max_iter_num + ) + +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() +# simulation.write_exact_solution_to_xdmf() +simulation.run(analyse_condition=analyse_condition) diff --git a/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/corner_subdomains/TP-TP-4-patch-pure-dd-corner_subdomains-avoiding-origin.py b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/corner_subdomains/TP-TP-4-patch-pure-dd-corner_subdomains-avoiding-origin.py new file mode 100755 index 0000000000000000000000000000000000000000..145014f69c24f4a2c1dc440b1fcf59d568ee047d --- /dev/null +++ b/TP-TP-2-patch-pure-dd-avoid-interface-at-origin/corner_subdomains/TP-TP-4-patch-pure-dd-corner_subdomains-avoiding-origin.py @@ -0,0 +1,570 @@ +#!/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 helpers as hlp +import datetime + +date = datetime.datetime.now() +datestr = date.strftime("%Y-%m-%d") +#import ufl as ufl + +# init sympy session +sym.init_printing() + +use_case = "TP-TP-4-patch-pure-dd-corner-subdomains-avoding-origin" +solver_tol = 5E-7 +max_iter_num = 1000 +FEM_Lagrange_degree = 1 + +############ GRID ####################### +mesh_resolution = 20 +timestep_size = 0.0001 +number_of_timesteps = 10000 +# 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 = 6 +starttime = 0 + +Lw = 0.25 #/timestep_size +Lnw=Lw + +lambda_w = 40 +lambda_nw = 40 + +include_gravity = False +debugflag = False +analyse_condition = True + +output_string = "./output/{}-{}_timesteps{}_polynomial_degree{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree) + +##### Domain and Interface #### +# global simulation domain domain +sub_domain0_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)] +# interface between subdomain1 and subdomain2 +interface12_vertices = [df.Point(0.6, 1.0), + df.Point(1.0, 0.6)] + +interface13_vertices = [df.Point(0.75, -1.0), + df.Point(1.0, -0.75)] + +interface14_vertices = [df.Point(-0.8, -1.0), + df.Point(-0.8, 1.0)] +# subdomain1. +sub_domain1_vertices = [interface14_vertices[0], + interface13_vertices[0], + interface13_vertices[1], + interface12_vertices[1], + interface12_vertices[0], + interface14_vertices[1] ] + +# 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: [sub_domain1_vertices[0], + sub_domain1_vertices[1]], + 1: [sub_domain1_vertices[2], + sub_domain1_vertices[3]], + 2: [sub_domain1_vertices[4], + sub_domain1_vertices[5]] +} +# subdomain2 +sub_domain2_vertices = [interface12_vertices[1], + sub_domain0_vertices[2], + interface12_vertices[0], + ] + +subdomain2_outer_boundary_verts = { + 0: sub_domain2_vertices +} + +# subdomain3 +sub_domain3_vertices = [interface13_vertices[0], + sub_domain0_vertices[1], + interface13_vertices[1], + ] + +subdomain3_outer_boundary_verts = { + 0: sub_domain3_vertices +} + +# subdomain4 +sub_domain4_vertices = [sub_domain0_vertices[0], + interface14_vertices[0], + interface14_vertices[1], + sub_domain0_vertices[3], + ] + +subdomain4_outer_boundary_verts = { + 0: [sub_domain4_vertices[2], + sub_domain4_vertices[3], + sub_domain4_vertices[0], + sub_domain4_vertices[1] + ] +} + +# 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, + sub_domain3_vertices, + sub_domain4_vertices + ] +# in the below list, index 0 corresponds to the 12 interface which has index 1 +interface_def_points = [interface12_vertices, + interface13_vertices, + interface14_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 +} + +# 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], + [1,3], + [1,4]] +isRichards = { + 1: False, # + 2: False, + 3: False, # + 4: False + } + + +viscosity = {# +# subdom_num : viscosity + 1: {'wetting' :1, + 'nonwetting': 1}, # + 2: {'wetting' :1, + 'nonwetting': 1}, + 3: {'wetting' :1, + 'nonwetting': 1}, # + 4: {'wetting' :1, + 'nonwetting': 1} +} + +porosity = {# +# subdom_num : porosity + 1: 1,# + 2: 1, + 3: 1, + 4: 1, +} + +# Dict of the form: { subdom_num : density } +densities = { + 1: {'wetting' :1, #997 + 'nonwetting': 1}, ##1225}, + 2: {'wetting' :1, + 'nonwetting': 1}, + 3: {'wetting' :1, + 'nonwetting': 1}, # + 4: {'wetting' :1, + 'nonwetting': 1} +} + +gravity_acceleration = 9.81 + + +L = {# +# subdom_num : subdomain L for L-scheme + 1 : {'wetting' :Lw, + 'nonwetting': Lnw},# + 2 : {'wetting' :Lw, + 'nonwetting': Lnw}, + 3 : {'wetting' :Lw, + 'nonwetting': Lnw},# + 4 : {'wetting' :Lw, + 'nonwetting': Lnw} +} + + +lambda_param = {# +# subdom_num : lambda parameter for the L-scheme + 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} +} + +## 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 + +_rel_perm1w = ft.partial(rel_perm1w) +_rel_perm1nw = ft.partial(rel_perm1nw) + +subdomain1_rel_perm = { + 'wetting': _rel_perm1w,# + 'nonwetting': _rel_perm1nw +} +## relative permeabilty functions on subdomain 2 +def rel_perm2w(s): + # relative permeabilty wetting on subdomain2 + return s**2 +def rel_perm2nw(s): + # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2 + return (1-s)**2 + +_rel_perm2w = ft.partial(rel_perm2w) +_rel_perm2nw = ft.partial(rel_perm2nw) + +subdomain2_rel_perm = { + 'wetting': _rel_perm2w,# + 'nonwetting': _rel_perm2nw +} +subdomain3_rel_perm = subdomain1_rel_perm +subdomain4_rel_perm = subdomain2_rel_perm +## dictionary of relative permeabilties on all domains. +relative_permeability = {# + 1: subdomain1_rel_perm, + 2: subdomain2_rel_perm, + 3: subdomain3_rel_perm, + 4: subdomain4_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 2*s + +def rel_perm2nw_prime(s): + # relative permeabilty on subdomain1 + return -2*(1-s) + +_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 +} + +subdomain3_rel_perm_prime = subdomain2_rel_perm_prime +subdomain4_rel_perm_prime = subdomain2_rel_perm_prime + +# dictionary of relative permeabilties on all domains. +ka_prime = { + 1: subdomain1_rel_perm_prime, + 2: subdomain2_rel_perm_prime, + 3: subdomain3_rel_perm_prime, + 4: subdomain4_rel_perm_prime, +} + + + +def saturation(pc, index): + # inverse capillary pressure-saturation-relationship + return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 1) + + + +def pc_sat_rel_sym(S, index): + # capillary pressure-saturation-relationship + return 1/S**(index+1) -1 + +pc_saturation_sym = { + 1: ft.partial(pc_sat_rel_sym, index=1), + 2: ft.partial(pc_sat_rel_sym, index=1), + 3: ft.partial(pc_sat_rel_sym, index=1), + 4: ft.partial(pc_sat_rel_sym, index=1), +} + +def saturation_sym(pc, index): + # inverse capillary pressure-saturation-relationship + return 1/((1 + pc)**(1/(index + 1))) + + +# derivative of S-pc relationship with respect to pc. This is needed for the +# construction of a analytic solution. +def saturation_sym_prime(pc, index): + # inverse capillary pressure-saturation-relationship + return -1/((index+1)*(1 + pc)**((index+2)/(index+1))) + + +# 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=1), + 2: ft.partial(saturation_sym, index=1), + 3: ft.partial(saturation_sym, index=1), + 4: ft.partial(saturation_sym, index=1) +} + +S_pc_sym_prime = { + 1: ft.partial(saturation_sym_prime, index=1), + 2: ft.partial(saturation_sym_prime, index=1), + 3: ft.partial(saturation_sym_prime, index=1), + 4: ft.partial(saturation_sym_prime, index=1) +} + +sat_pressure_relationship = { + 1: ft.partial(saturation, index=1), + 2: ft.partial(saturation, index=1), + 3: ft.partial(saturation, index=1), + 4: ft.partial(saturation, index=1) +} + + +############################################# +# Manufacture source expressions with sympy # +############################################# +x, y = sym.symbols('x[0], x[1]') # needed by UFL +t = sym.symbols('t', positive=True) + +symbols = { "x": x, + "y": y, + "t": t} + +# epsilon_x_inner = 0.7 +# epsilon_x_outer = 0.99 +# epsilon_y_inner = epsilon_x_inner +# epsilon_y_outer = epsilon_x_outer +# +# def mollifier(x, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x/epsilon)**2) + 1) +# return out_expr +# +# mollifier_handle = ft.partial(mollifier, epsilon=epsilon_x_inner) +# +# pw_sym_x = sym.Piecewise( +# (mollifier_handle(x), x**2 < epsilon_x_outer**2), +# (0, True) +# ) +# pw_sym_y = sym.Piecewise( +# (mollifier_handle(y), y**2 < epsilon_y_outer**2), +# (0, True) +# ) +# +# def mollifier2d(x, y, epsilon): +# """ one d mollifier """ +# out_expr = sym.exp(-1/(1-(x**2 + y**2)/epsilon**2) + 1) +# return out_expr +# +# mollifier2d_handle = ft.partial(mollifier2d, epsilon=epsilon_x_outer) +# +# pw_sym2d_x = sym.Piecewise( +# (mollifier2d_handle(x, y), x**2 + y**2 < epsilon_x_outer**2), +# (0, True) +# ) +# +# zero_on_epsilon_shrinking_of_subdomain = sym.Piecewise( +# (mollifier_handle(sym.sqrt(x**2 + y**2)+2*epsilon_x_inner), ((-2*epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<=epsilon_x_inner))), +# (mollifier_handle(sym.sqrt(x**2 + y**2)-2*epsilon_x_inner), ((epsilon_x_inner<sym.sqrt(x**2 + y**2)) & (sym.sqrt(x**2 + y**2)<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_x = sym.Piecewise( +# (mollifier_handle(x+2*epsilon_x_inner), ((-2*epsilon_x_inner<x) & (x<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=x) & (x<=epsilon_x_inner))), +# (mollifier_handle(x-2*epsilon_x_inner), ((epsilon_x_inner<x) & (x<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_epsilon_shrinking_of_subdomain_y = sym.Piecewise( +# (1, y<=-2*epsilon_x_inner), +# (mollifier_handle(y+2*epsilon_x_inner), ((-2*epsilon_x_inner<y) & (y<-epsilon_x_inner))), +# (0, ((-epsilon_x_inner<=y) & (y<=epsilon_x_inner))), +# (mollifier_handle(y-2*epsilon_x_inner), ((epsilon_x_inner<y) & (y<2*epsilon_x_inner))), +# (1, True), +# ) +# +# zero_on_shrinking = zero_on_epsilon_shrinking_of_subdomain #zero_on_epsilon_shrinking_of_subdomain_x + zero_on_epsilon_shrinking_of_subdomain_y +# gaussian = pw_sym2d_x# pw_sym_y*pw_sym_x +# cutoff = gaussian/(gaussian + zero_on_shrinking) +# +# +# sat_sym = { +# 1: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t), +# 2: 0.5 + 0.25*sym.sin(x-t)*sym.cos(y-t) +# } +# +# Spc = { +# 1: sym.Piecewise((pc_saturation_sym[1](sat_sym[1]), sat_sym[1] > 0), (pc_saturation_sym[1](sat_sym[1]), 1>=sat_sym[1]), (0, True)), +# 2: sym.Piecewise((pc_saturation_sym[2](sat_sym[2]), sat_sym[2] > 0), (pc_saturation_sym[2](sat_sym[2]), 2>=sat_sym[2]), (0, True)) +# } +# +# p1w = (-1 - (1+t*t)*(1 + x*x + y*y))#*cutoff +# p2w = p1w +# p_e_sym = { +# 1: {'wetting': p1w, +# 'nonwetting': (p1w + Spc[1])}, #*cutoff}, +# 2: {'wetting': p2w, +# 'nonwetting': (p2w + Spc[2])}, #*cutoff}, +# } + +p_e_sym = { + 1: {'wetting': (-6 - (1+t*t)*(1 + x*x + y**2)), #*cutoff, + 'nonwetting': -1 -t*(1.1 + x**2)}, #*(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**2)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2, + 'nonwetting': -1 -t*(1.1 + x**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, + 3: {'wetting': (-6 - (1+t*t)*(1 + x*x + y**2)), #*cutoff, + 'nonwetting': -1 -t*(1.1 + x**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, + 4: {'wetting': (-6 - (1+t*t)*(1 + x*x + y**2)), #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2, + 'nonwetting': -1 -t*(1.1 + x**2)}, #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2}, +} + + +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']}) + + + +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]} + ) + + +# def saturation(pressure, subdomain_index): +# # inverse capillary pressure-saturation-relationship +# return df.conditional(pressure < 0, 1/((1 - pressure)**(1/(subdomain_index + 1))), 1) +# +# sa + +write_to_file = { + 'meshes_and_markers': True, + 'L_iterations': True +} + + +# initialise LDD simulation class +simulation = ldd.LDDsimulation( + tol=1E-14, + LDDsolver_tol=solver_tol, + debug=debugflag, + max_iter_num=max_iter_num, + FEM_Lagrange_degree=FEM_Lagrange_degree + ) + +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() +# simulation.write_exact_solution_to_xdmf() +simulation.run(analyse_condition=analyse_condition)