TP-multi-patch-with-gravity-same-wetting-phase-as-RR.py

#!/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
import os
import pandas as pd

date = datetime.datetime.now()
datestr = date.strftime("%Y-%m-%d")

# init sympy session
sym.init_printing()
# solver_tol = 6E-7
use_case = "TP-TP-4-patch-nonwetting-zero"
max_iter_num = 1000
FEM_Lagrange_degree = 1
mesh_study = False
resolutions = {
                # 1: 1e-7,  # h=2
                # 2: 2e-5,  # h=1.1180
                # 4: 1e-6,  # h=0.5590
                # 8: 1e-6,  # h=0.2814
                # 16: 5e-7, # h=0.1412
                32: 7e-7,   # h=0.0706
                # 64: 1e-6,
                # 128: 5e-7
                }

############ GRID #######################
# mesh_resolution = 20
timestep_size = 0.005
number_of_timesteps = 250
plot_timestep_every = 1
# 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 = 2
starttime = 0.0

Lw = 0.05 #/timestep_size
Lnw=Lw

lambda_w = 20
lambda_nw = 20

include_gravity = True
debugflag = False
analyse_condition = True

if mesh_study:
    output_string = "./output/{}-{}_timesteps{}_P{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree)
else:
    for tol in resolutions.values():
        solver_tol = tol
    output_string = "./output/{}-{}_timesteps{}_P{}_solver_tol{}".format(datestr, use_case, number_of_timesteps, FEM_Lagrange_degree, solver_tol)

# toggle what should be written to files
if mesh_study:
    write_to_file = {
        'space_errornorms': True,
        'meshes_and_markers': True,
        'L_iterations_per_timestep': False,
        'solutions': False,
        'absolute_differences': False,
        'condition_numbers': analyse_condition,
        'subsequent_errors': False
    }
else:
    write_to_file = {
        'space_errornorms': True,
        'meshes_and_markers': True,
        'L_iterations_per_timestep': False,
        'solutions': True,
        'absolute_differences': True,
        'condition_numbers': analyse_condition,
        'subsequent_errors': True
    }
# ----------------------------------------------------------------------------#
# ------------------- 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)]
# interfaces
interface12_vertices = [df.Point(0.0, 0.0),
                        df.Point(1.0, 0.0)]

interface14_vertices = [df.Point(0.0, 0.0),
                        df.Point(0.0, 1.0)]

interface23_vertices = [df.Point(0.0, 0.0),
                        df.Point(0.0, -1.0)]

interface34_vertices = [df.Point(-1.0, 0.0),
                        df.Point(0.0, 0.0)]
# subdomain1.
sub_domain1_vertices = [interface12_vertices[0],
                        interface12_vertices[1],
                        sub_domain0_vertices[2],
                        df.Point(0.0, 1.0)]

# 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 inter-
# nal, the dictionary entry should be 0: None
subdomain1_outer_boundary_verts = {
    0: [interface12_vertices[1],
        sub_domain0_vertices[2],
        df.Point(0.0, 1.0)]
}
# subdomain2
sub_domain2_vertices = [interface23_vertices[1],
                        sub_domain0_vertices[1],
                        interface12_vertices[1],
                        interface12_vertices[0]]

subdomain2_outer_boundary_verts = {
    0: [df.Point(0.0, -1.0),
        sub_domain0_vertices[1],
        interface12_vertices[1]]
}
sub_domain3_vertices = [interface34_vertices[0],
                        sub_domain0_vertices[0],
                        interface23_vertices[1],
                        interface23_vertices[0]]

subdomain3_outer_boundary_verts = {
    0: [interface34_vertices[0],
        sub_domain0_vertices[0],
        interface23_vertices[1]]
}

sub_domain4_vertices = [interface34_vertices[0],
                        interface34_vertices[1],
                        interface14_vertices[1],
                        sub_domain0_vertices[3]]

subdomain4_outer_boundary_verts = {
    0: [interface14_vertices[1],
        sub_domain0_vertices[3],
        interface34_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,
                        sub_domain3_vertices,
                        sub_domain4_vertices]
# in the below list, index 0 corresponds to the 12 interface which has global
# marker value 1
interface_def_points = [interface12_vertices,
                        interface14_vertices,
                        interface23_vertices,
                        interface34_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], [1, 4], [2, 3], [3, 4]]

# 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
}

isRichards = {
    1: False,
    2: False,
    3: False,
    4: False
    }

viscosity = {
    1: {'wetting' :1,
         'nonwetting': 1},
    2: {'wetting' :1,
         'nonwetting': 1},
    3: {'wetting' :1,
         'nonwetting': 1},
    4: {'wetting' :1,
         'nonwetting': 1},
}

# Dict of the form: { subdom_num : density }
densities = {
    1: {'wetting': 1,  #997
         'nonwetting':1},  #1.225}},
    2: {'wetting': 1,  #997
         'nonwetting':1},  #1.225}},
    3: {'wetting': 1,  #997
         'nonwetting':1},  #1.225}},
    4: {'wetting': 1,  #997
         'nonwetting':1},  #1.225}}
}

gravity_acceleration = 9.81
# porosities taken from
# https://www.geotechdata.info/parameter/soil-porosity.html
# Dict of the form: { subdom_num : porosity }
porosity = {
    1: 1,  #0.2,  # Clayey gravels, clayey sandy gravels
    2: 1,  #0.22, # Silty gravels, silty sandy gravels
    3: 1,  #0.37, # Clayey sands
    4: 1,  #0.2 # Silty or sandy clay
}

# subdom_num : subdomain L for L-scheme
L = {
    1: {'wetting' :Lw,
         'nonwetting': Lnw},
    2: {'wetting' :Lw,
         'nonwetting': Lnw},
    3: {'wetting' :Lw,
         'nonwetting': Lnw},
    4: {'wetting' :Lw,
         'nonwetting': Lnw}
}

# subdom_num : lambda parameter for the L-scheme
lambda_param = {
    1: {'wetting': lambda_w,
         'nonwetting': lambda_nw},#
    2: {'wetting': lambda_w,
         'nonwetting': lambda_nw},#
    3: {'wetting': lambda_w,
         'nonwetting': lambda_nw},#
    4: {'wetting': lambda_w,
         'nonwetting': lambda_nw},#
}


# relative permeabilty functions on subdomain 1
def rel_perm1w(s):
    # relative permeabilty on subdomain1
    return s**2


def rel_perm1nw(s):
    # relative permeabilty nonwetting on subdomain1
    return (1-s)**2


## relative permeabilty functions on subdomain 2
# relative permeabilty functions on subdomain 2
def rel_perm2w(s):
    # relative permeabilty on subdomain2
    return s**3


def rel_perm2nw(s):
    # relative permeabilty nonwetting on subdosym.cos(0.8*t - (0.8*x + 1/7*y))main2
    return (1-s)**3


_rel_perm1w = ft.partial(rel_perm1w)
_rel_perm1nw = ft.partial(rel_perm1nw)
_rel_perm2w = ft.partial(rel_perm2w)
_rel_perm2nw = ft.partial(rel_perm2nw)

subdomain1_rel_perm = {
    'wetting': _rel_perm1w,#
    'nonwetting': _rel_perm1nw
}

subdomain2_rel_perm = {
    'wetting': _rel_perm2w,#
    'nonwetting': _rel_perm2nw
}


subdomain3_rel_perm = subdomain2_rel_perm.copy()
subdomain4_rel_perm = subdomain1_rel_perm.copy()

# 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 3*s**2

def rel_perm2nw_prime(s):
    # relative permeabilty on subdomain1
    return -3*(1-s)**2

_rel_perm1w_prime = ft.partial(rel_perm1w_prime)
_rel_perm1nw_prime = ft.partial(rel_perm1nw_prime)
_rel_perm2w_prime = ft.partial(rel_perm2w_prime)
_rel_perm2nw_prime = ft.partial(rel_perm2nw_prime)

subdomain1_rel_perm_prime = {
    'wetting': _rel_perm1w_prime,
    'nonwetting': _rel_perm1nw_prime
}


subdomain2_rel_perm_prime = {
    'wetting': _rel_perm2w_prime,
    'nonwetting': _rel_perm2nw_prime
}

# _rel_perm3_prime = ft.partial(rel_perm2_prime)
subdomain3_rel_perm_prime = subdomain2_rel_perm_prime.copy()

# _rel_perm4_prime = ft.partial(rel_perm1_prime)
subdomain4_rel_perm_prime = subdomain1_rel_perm_prime.copy()

# 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
}


# this function needs to be monotonically decreasing in the capillary_pressure.
# since in the richards case pc=-pw, this becomes as a function of pw a mono
# tonically INCREASING function like in our Richards-Richards paper. However
# since we unify the treatment in the code for Richards and two-phase, we need
# the same requierment
# for both cases, two-phase and Richards.
def saturation(pc, index):
    # inverse capillary pressure-saturation-relationship
    return df.conditional(pc > 0, 1/((1 + pc)**(1/(index + 1))), 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=2),
    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=2),
    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=2),
    3: ft.partial(saturation, index=2),
    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)

p_e_sym = {
    1: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
        'nonwetting': 0.0*t},
    2: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
        'nonwetting': 0.0*t},
    3: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
        'nonwetting': 0.0*t},
    4: {'wetting': (-3.0 - (1.0 + t*t)*(1.0 + x*x + y*y)),  #*(sym.sin((1+y)/2*sym.pi)*sym.sin((1+x)/2*sym.pi))**2,
        'nonwetting': 0.0*t}
}

# pc_e_sym = {
#     1: -1*p_e_sym[1]['wetting'],
#     2: -1*p_e_sym[2]['wetting'],
#     3: -1*p_e_sym[3]['wetting'],
#     4: -1*p_e_sym[4]['wetting']
# }

pc_e_sym = dict()
for subdomain, isR in isRichards.items():
    if isR:
        pc_e_sym.update({subdomain: -p_e_sym[subdomain]['wetting']})
    else:
        pc_e_sym.update({subdomain: p_e_sym[subdomain]['nonwetting']
                                        - p_e_sym[subdomain]['wetting']})

symbols = {"x": x,
           "y": y,
           "t": t}
# turn above symbolic code into exact solution for dolphin and
# construct the rhs that matches the above exact solution.
exact_solution_example = hlp.generate_exact_solution_expressions(
                        symbols=symbols,
                        isRichards=isRichards,
                        symbolic_pressure=p_e_sym,
                        symbolic_capillary_pressure=pc_e_sym,
                        saturation_pressure_relationship=S_pc_sym,
                        saturation_pressure_relationship_prime=S_pc_sym_prime,
                        viscosity=viscosity,
                        porosity=porosity,
                        relative_permeability=relative_permeability,
                        relative_permeability_prime=ka_prime,
                        densities=densities,
                        gravity_acceleration=gravity_acceleration,
                        include_gravity=include_gravity,
                        )
source_expression = exact_solution_example['source']
exact_solution = exact_solution_example['exact_solution']
initial_condition = exact_solution_example['initial_condition']


# Dictionary of dirichlet boundary conditions.
dirichletBC = dict()
# similarly to the outer boundary dictionary, if a patch has no outer boundary
# None should be written instead of an expression.
# This is a bit of a brainfuck:
# dirichletBC[ind] gives a dictionary of the outer boundaries of subdomain ind.
# Since a domain patch can have several disjoint outer boundary parts, the
# expressions need to get an enumaration index which starts at 0.
# So dirichletBC[ind][j] is the dictionary of outer dirichlet conditions of
# subdomain ind and boundary part j.
# Finally, dirichletBC[ind][j]['wetting'] and dirichletBC[ind][j]['nonwetting']
# return the actual expression needed for the dirichlet condition for both
# phases if present.

# subdomain index: {outer boudary part index: {phase: expression}}
for subdomain in isRichards.keys():
    # if subdomain has no outer boundary, outer_boundary_def_points[subdomain] is None
    if outer_boundary_def_points[subdomain] is None:
        dirichletBC.update({subdomain: None})
    else:
        dirichletBC.update({subdomain: dict()})
        # set the dirichlet conditions to be the same code as exact solution on
        # the subdomain.
        for outer_boundary_ind in outer_boundary_def_points[subdomain].keys():
            dirichletBC[subdomain].update(
                {outer_boundary_ind: exact_solution[subdomain]}
                )



for mesh_resolution, solver_tol in resolutions.items():
    # 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,
        mesh_study=mesh_study
        )

    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,
                              plot_timestep_every=plot_timestep_every,
                              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()
    output_dir = simulation.output_dir
    # simulation.write_exact_solution_to_xdmf()
    output = simulation.run(analyse_condition=analyse_condition)
    for subdomain_index, subdomain_output in output.items():
        mesh_h = subdomain_output['mesh_size']
        for phase, different_errornorms in subdomain_output['errornorm'].items():
            filename = output_dir + "subdomain{}-space-time-errornorm-{}-phase.csv".format(subdomain_index, phase)
            # for errortype, errornorm in different_errornorms.items():

                # eocfile = open("eoc_filename", "a")
                # eocfile.write( str(mesh_h) + " " + str(errornorm) + "\n" )
                # eocfile.close()
                # if subdomain.isRichards:mesh_h
            data_dict = {
                'mesh_parameter': mesh_resolution,
                'mesh_h': mesh_h,
            }
            for error_type, errornorms in different_errornorms.items():
                data_dict.update(
                    {error_type: errornorms}
                )
            errors = pd.DataFrame(data_dict, index=[mesh_resolution])
            # check if file exists
            if os.path.isfile(filename) == True:
                with open(filename, 'a') as f:
                    errors.to_csv(f, header=False, sep='\t', encoding='utf-8', index=False)
            else:
                errors.to_csv(filename, sep='\t', encoding='utf-8', index=False)