diff --git a/LDDsimulation/LDDsimulation.py b/LDDsimulation/LDDsimulation.py
index 0707f040812f5a311e5a2ef2104207c0e2fb4a12..85ec4c138c0d94a851d0124fbc7867aed1f18aa3 100644
--- a/LDDsimulation/LDDsimulation.py
+++ b/LDDsimulation/LDDsimulation.py
@@ -11,6 +11,7 @@ import domainPatch as dp
import helpers as hlp
import pandas as pd
import sys
+import copy
# import errno
import h5py
from termcolor import colored
@@ -180,7 +181,7 @@ class LDDsimulation(object):
starttime: float,#
timestep_size: tp.Dict[int, float],#
number_of_timesteps: int,#
- number_of_timesteps_to_analise: int,#
+ number_of_timesteps_to_analyse: int,#
sources: tp.Dict[int, tp.Dict[str, str]],#
initial_conditions: tp.Dict[int, tp.Dict[str, str]],#
dirichletBC_Expressions: tp.Dict[int, tp.Dict[str, str]],#
@@ -207,7 +208,7 @@ class LDDsimulation(object):
self.t = starttime
# total number of timesteps.
self.number_of_timesteps = number_of_timesteps
- self.number_of_timesteps_to_analise = number_of_timesteps_to_analise
+ self.number_of_timesteps_to_analyse = number_of_timesteps_to_analyse
self.timestep_size = timestep_size
self.sources = sources
self.initial_conditions = initial_conditions
@@ -239,6 +240,9 @@ class LDDsimulation(object):
# initialise exact solution expression dictionary if we have an exact
# solution case.
self._init_exact_solution_expression()
+ if self.exact_solution:
+ print("writing exact solution for all time steps to xdmf files ...")
+ self.write_exact_solution_to_xdmf()
self._initialised = True
def run(self, analyse_solver_at_timesteps: tp.List[float] = None):
@@ -589,7 +593,7 @@ class LDDsimulation(object):
Lam = subdomain.lambda_param['wetting']
name1 ="subdomain" + "{}".format(subdom_ind)
# only show three digits of the mesh size.
- name2 = "_dt" + "{}".format(tau)+"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)
+ name2 = "_dt" + "{}".format(tau)+"_hmax" + "{number:.{digits}f}".format(number=h, digits=3)
name3 = "_lambda" + "{}".format(Lam) + "_Lw" + "{}".format(L["wetting"])
filename_param_part = name1 + name2 + name3
filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
@@ -598,12 +602,12 @@ class LDDsimulation(object):
Lam_nw = subdomain.lambda_param['nonwetting']
filename_param_part = "subdomain" + "{}".format(subdom_ind)\
+"_dt" + "{}".format(tau)\
- +"_hmin" + "{number:.{digits}f}".format(number=h, digits=3)\
+ +"_hmax" + "{number:.{digits}f}".format(number=h, digits=3)\
+"_lambda_w" + "{}".format(Lam_w)\
+"_lambda_nw" + "{}".format(Lam_nw)\
+"_Lw" + "{}".format(L["wetting"])\
+"_Lnw" + "{}".format(L["nonwetting"])
- filename = self.output_dir+"_solution" +".xdmf"
+ filename = self.output_dir + filename_param_part + "_solution" +".xdmf"
# create solution files and populate dictionary.
self.output_filename_parameter_part.update(
{subdom_ind: filename_param_part}
@@ -612,6 +616,22 @@ class LDDsimulation(object):
{subdom_ind: SolutionFile(self._mpi_communicator, filename)}
)
+ def write_exact_solution_to_xdmf(self):
+ """
+ evaluate the exact solution of the simulation (in case there is one) at
+ all timesteps and write it to the solution file.
+ """
+ t = copy.copy(self.starttime)
+ for timestep in range(0,self.number_of_timesteps+1):
+ for subdom_ind, subdomain in self.subdomain.items():
+ file = self.solution_file[subdom_ind]
+ for phase in subdomain.has_phases:
+ pa_exact = subdomain.pressure_exact[phase]
+ pa_exact.t = t
+ tmp_exact_pressure = df.interpolate(pa_exact, subdomain.function_space[phase])
+ tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
+ file.write(tmp_exact_pressure, t)
+ t += self.timestep_size
def write_subsequent_errors_to_csv(self,#
filename: str, #
@@ -1074,24 +1094,24 @@ class LDDsimulation(object):
by setting
"""
# We want to write out csv files containing subsequent L2 errors of iterations for
- # number_of_timesteps_to_analise many timesteps. We want to do this with
- # analyse_timesteps = np.linspace(starttime, number_of_timesteps, number_of_timesteps_to_analise)
- # see below. Therefor, we need to check if (number_of_timesteps_to_analise-1) is
+ # number_of_timesteps_to_analyse many timesteps. We want to do this with
+ # analyse_timesteps = np.linspace(starttime, number_of_timesteps, number_of_timesteps_to_analyse)
+ # see below. Therefor, we need to check if (number_of_timesteps_to_analyse-1) is
# a divisor of number_of_timesteps,
- if self.number_of_timesteps_to_analise == 0:
+ if self.number_of_timesteps_to_analyse == 0:
self.analyse_timesteps = None
- elif self.number_of_timesteps_to_analise == 1:
+ elif self.number_of_timesteps_to_analyse == 1:
self.analyse_timesteps = [0]
else:
- if not self.number_of_timesteps%(self.number_of_timesteps_to_analise-1) == 0:
- # if number_of_timesteps%(number_of_timesteps_to_analise-1) is not 0, we
+ if not self.number_of_timesteps%(self.number_of_timesteps_to_analyse-1) == 0:
+ # if number_of_timesteps%(number_of_timesteps_to_analyse-1) is not 0, we
# round the division error to the nearest integer and redefine number_of_timesteps
- print(f"since number_of_timesteps/number_of_timesteps_to_analise = {self.number_of_timesteps/self.number_of_timesteps_to_analise}",
+ print(f"since number_of_timesteps/number_of_timesteps_to_analyse = {self.number_of_timesteps/self.number_of_timesteps_to_analyse}",
"is not an integer,\n")
- new_analysing_interval = int(round(self.number_of_timesteps/(self.number_of_timesteps_to_analise-1)))
- self.number_of_timesteps = new_analysing_interval*(self.number_of_timesteps_to_analise)
+ new_analysing_interval = int(round(self.number_of_timesteps/(self.number_of_timesteps_to_analyse-1)))
+ self.number_of_timesteps = new_analysing_interval*(self.number_of_timesteps_to_analyse)
print(f"I'm resetting number_of_timesteps = {self.number_of_timesteps}\n")
- self.analyse_timesteps = np.linspace(0, self.number_of_timesteps, self.number_of_timesteps_to_analise, dtype=int)
+ self.analyse_timesteps = np.linspace(0, self.number_of_timesteps, self.number_of_timesteps_to_analyse, dtype=int)
print(f"the following timesteps will be analysed: {self.analyse_timesteps}")
self.endtime = self.number_of_timesteps*self.timestep_size
diff --git a/LDDsimulation/domainPatch.py b/LDDsimulation/domainPatch.py
index 143780928fa999d0cc9d935fd14223055fdb4589..ad5411327c16e7ae83e4ae230a714d62d1d875d6 100644
--- a/LDDsimulation/domainPatch.py
+++ b/LDDsimulation/domainPatch.py
@@ -20,7 +20,9 @@ import numpy as np
# import domainPatch as dp
import typing as tp
import ufl as fl
+import LDDsimulation
# Interface = tp.NewType('Interface', tp.any)
+# SolutionFile = tp.NewType('LDDsimulation.SolutionFile', tp.object)
#import domainPatch as domainPatch
@@ -785,9 +787,10 @@ class DomainPatch(df.SubDomain):
)
- def write_solution_to_xdmf(self, file: str, #
+ def write_solution_to_xdmf(self, file: tp.Type[LDDsimulation.SolutionFile], #
time: float = None,#
write_iter_for_fixed_time: bool = False,#
+ exact_solution: bool = False,#
):
"""
Weither write the solution of the simulation and corresponding time
@@ -811,6 +814,11 @@ class DomainPatch(df.SubDomain):
if not write_iter_for_fixed_time:
self.pressure[phase].rename("pressure_"+"{}".format(phase), "pressure_"+"{}".format(phase))
file.write(self.pressure[phase], t)
+ # if self.pressure_exact:
+ # self.pressure_exact[phase].t = t
+ # tmp_exact_pressure = df.interpolate(self.pressure_exact[phase], self.function_space[phase])
+ # tmp_exact_pressure.rename("exact_pressure_"+"{}".format(phase), "exactpressure_"+"{}".format(phase))
+ # file.write(tmp_exact_pressure, t)
else:
i = self.iteration_number
self.pressure[phase].rename("pressure_"+"{}".format(phase), #
@@ -828,28 +836,6 @@ class DomainPatch(df.SubDomain):
self.gli_function[phase].rename("gli_function_"+"{}".format(phase),#
"all_gli_terms(t="+"{}".format(t)+")_"+"{}".format(phase))
file.write(self.gli_function[phase], i)
- # # write unknowns to file
- # rho.rename("rho","density")
- # v.rename("v","velocity")
- # phi.rename("phi","phasefield")
- # mu.rename("mu","dfdc")
- # tau.rename("tau","dfdrho")
- # sigma.rename("sigma","gradc")
- #
- # for var in [rho, v, phi ,mu, tau, sigma]:
- # file.write(var, self.t)
- #
- # # write pressure to file
- # p = df.project(self.pressure(),df.FunctionSpace(self._mesh, self._element.sub_elements()[0]),solver_type="bicgstab")
- # if self.project_to_cg == True:
- # p = df.project(self.pressure(),df.FunctionSpace(self._mesh, self._element.sub_elements()[0].reconstruct(family="CG")),solver_type="bicgstab")
- #
- # p.rename("p","pressure")
- # file.write(p, self.t)
-
-
-
-
# END OF CLASS DomainPatch
diff --git a/TP-TP-patch-test-case/TP-TP-2-patch-test.py b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
index dff4770c246d1afdcc89d561b20d144b57034019..7ad8c9719a3d987a3a8f57bccc0662b39df92671 100755
--- a/TP-TP-patch-test-case/TP-TP-2-patch-test.py
+++ b/TP-TP-patch-test-case/TP-TP-2-patch-test.py
@@ -89,12 +89,12 @@ isRichards = {
############ GRID ########################ΓΌ
-mesh_resolution = 20
-timestep_size = 4*0.001
-number_of_timesteps = 500
+mesh_resolution = 100
+timestep_size = 2*0.001
+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_analise = 11
+number_of_timesteps_to_analyse = 11
starttime = 0
viscosity = {#
@@ -122,9 +122,9 @@ L = {#
lambda_param = {#
# subdom_num : lambda parameter for the L-scheme
1 : {'wetting' :140,
- 'nonwetting': 2500},#
+ 'nonwetting': 1000},#
2 : {'wetting' :140,
- 'nonwetting': 2500}
+ 'nonwetting': 1000}
}
## relative permeabilty functions on subdomain 1
@@ -334,7 +334,7 @@ simulation.set_parameters(output_dir = "./output/",#
saturation = S_pc_rel,#
starttime = starttime,#
number_of_timesteps = number_of_timesteps,
- number_of_timesteps_to_analise = number_of_timesteps_to_analise,
+ number_of_timesteps_to_analyse = number_of_timesteps_to_analyse,
timestep_size = timestep_size,#
sources = source_expression,#
initial_conditions = initial_condition,#
@@ -344,4 +344,5 @@ simulation.set_parameters(output_dir = "./output/",#
)
simulation.initialise()
+# simulation.write_exact_solution_to_xdmf()
simulation.run()