1005 lines
48 KiB
Python
1005 lines
48 KiB
Python
"""
|
|
This module contains the problem definition and solver for iNSE ALE article.
|
|
It enables to replicate all results
|
|
"""
|
|
__version__ = "0.2"
|
|
__author__ = "Reidmen Arostica"
|
|
__date__ = "18/07/2020"
|
|
|
|
# Libraries
|
|
import fenics as fn
|
|
from dolfin import __version__ as df_version
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
import time
|
|
from common import *
|
|
|
|
# Basic info to user
|
|
print("DOLFIN VERSION: {}".format(df_version))
|
|
print("SOLVER VERSION: {}".format(__version__))
|
|
print("__author__: {}".format(__author__))
|
|
|
|
|
|
def iNSE_Monolithic_Solver(parameters, deformation, stabilization, solver_to_use='LU'):
|
|
""" This solver computes the incompressible Navier-Stokes (iNSE) problem in the ALE formalism.
|
|
To do so, the FEniCS framework is used and the convective term within iNSE is linealized.
|
|
To solve the problem, a monolithic scheme is employed with finite diferences for time
|
|
and FE for space.
|
|
...
|
|
|
|
Input Parameters
|
|
----------------
|
|
parameter: dict
|
|
dictionary containing several options for the scheme used.
|
|
deformation: dict
|
|
dictionary containing deformation mapping properties and options.
|
|
stabilization: dict
|
|
dictionary specifying stabilization used.
|
|
|
|
|
|
Output
|
|
----------------
|
|
None.
|
|
The solver print the results to console, and returns None, saving solutions and figures.
|
|
"""
|
|
# Starting time
|
|
start = time.time()
|
|
# Default parameters
|
|
mu, rho_f = parameters['mu'], parameters['rho_f']
|
|
Nx, Ny = parameters['Nx'], parameters['Ny']
|
|
tau, T = parameters['tau'], parameters['T']
|
|
|
|
def def_mesh_tuple(def_params: dict) -> tuple:
|
|
"""
|
|
Computes a list (len 2) containing the parametrized mesh velocity.
|
|
...
|
|
|
|
Input Parameters
|
|
----------------
|
|
def_params: dict
|
|
dictionary with keys 'expression' and 'output_name'.
|
|
|
|
"""
|
|
if 'expression' in def_params.keys():
|
|
if 'output_name' in def_params.keys():
|
|
output_name = def_params['output_name']
|
|
else:
|
|
output_name = ('X', 'Y')
|
|
return def_params['expression'], output_name
|
|
else:
|
|
raise Exception("Key 'expression' was not found in the dictionary!")
|
|
|
|
def exps_at_bdry(bdry_params: dict, def_params):
|
|
"""
|
|
Computes the boundary FEniCS expressions for this particular Dirichlet Homogeneous case.
|
|
"""
|
|
u_inlet = fn.Constant((0.0, 0.0))
|
|
u_outlet = fn.Constant((0.0, 0.0))
|
|
u_bottom = fn.Constant((0.0, 0.0))
|
|
u_top = fn.Constant((0.0, 0.0))
|
|
return u_inlet, u_outlet, u_bottom, u_top
|
|
|
|
# Define boundary domain classes
|
|
class inlet(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[0], 0) and on_boundary
|
|
class outlet(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[0], 6) and on_boundary
|
|
class top(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[1], 1) and on_boundary
|
|
class bottom(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[1], -1) and on_boundary
|
|
|
|
# Define mesh object
|
|
mesh_f = fn.RectangleMesh(fn.Point(0, -1), fn.Point(6, 1), Nx, Ny)
|
|
#h = fn.CellDiameter(mesh_f)
|
|
|
|
# Obtain mesh deformation tuple and output name
|
|
def_tuple, def_name = def_mesh_tuple(deformation)
|
|
|
|
# Define displacement expressions at previous and actual time stepss
|
|
def_n = fn.Expression(def_tuple, t=0.0, domain=mesh_f, degree=2)
|
|
def_t = fn.Expression(def_tuple, t=tau, domain=mesh_f, degree=2)
|
|
|
|
# Functional Spaces (Mixed Element)
|
|
U = fn.VectorElement("CG", mesh_f.ufl_cell(), parameters['deg_U'])
|
|
P = fn.FiniteElement("CG", mesh_f.ufl_cell(), parameters['deg_P'])
|
|
M = fn.FunctionSpace(mesh_f, fn.MixedElement([U, P]))
|
|
|
|
# Functional spaces (Separately)
|
|
U_ = fn.VectorFunctionSpace(mesh_f, "CG", parameters['deg_U'])
|
|
P_ = fn.FunctionSpace(mesh_f, "CG", parameters['deg_P'])
|
|
|
|
# Assigner between spaces
|
|
assigner_to_M = fn.FunctionAssigner(M, [U_, P_])
|
|
assigner_to_UP = fn.FunctionAssigner([U_, P_], M)
|
|
|
|
# Boundary expressions in the reference domain
|
|
u_inflow, u_outflow, u_bottom, u_top = exps_at_bdry(parameters, deformation)
|
|
# Define and mark boundaries
|
|
bdry_markers = fn.MeshFunction('size_t', mesh_f, mesh_f.topology().dim()-1)
|
|
inlet().mark(bdry_markers, 1) # Inlet marked with label 1
|
|
outlet().mark(bdry_markers, 2) # Outlet marked with label 2
|
|
bottom().mark(bdry_markers, 3) # Bottom marked with label 3
|
|
top().mark(bdry_markers, 4) # Top marked with label 4
|
|
|
|
# Define measures for weak formulation
|
|
dx_f = fn.Measure('dx', domain=mesh_f, subdomain_data=bdry_markers)
|
|
#ds_f = fn.Measure('ds', domain=mesh_f, subdomain_data=bdry_markers)
|
|
|
|
# Define boundary conditions with pressure condition to obtain uniqueness (case LU solver!)
|
|
bdries = [fn.DirichletBC(M.sub(0), u_inflow, bdry_markers, 1),
|
|
fn.DirichletBC(M.sub(0), u_outflow, bdry_markers, 2),
|
|
fn.DirichletBC(M.sub(0), u_bottom, bdry_markers, 3),
|
|
fn.DirichletBC(M.sub(0), u_top, bdry_markers, 4),
|
|
fn.DirichletBC(M.sub(1), fn.Constant(0.), "near(x[0], 0) && near(x[1], 0)", 'pointwise')]
|
|
|
|
# Trial/Test function definition
|
|
(u, p) = fn.TrialFunctions(M)
|
|
(v, q) = fn.TestFunctions(M)
|
|
aux_n = fn.Function(M, name='MixedSolution')
|
|
# Define initial condition (non-zero)
|
|
if parameters["init_cond"] == True:
|
|
print("Using a non-zero initial condition (2D-Gaussian) for the velocity field!")
|
|
cpp_aux = '0.001*(1 - pow(x[1], 2))*x[0]*(6 - x[0])'
|
|
aux_n.assign(fn.Expression((cpp_aux, '0.', '0.'), domain=mesh_f, degree=2))
|
|
u_n, _ = fn.split(aux_n)
|
|
else:
|
|
u_n, _ = fn.split(aux_n)
|
|
|
|
# Define normal for backflow stabilization term
|
|
#N = fn.FacetNormal(mesh_f)
|
|
# Compute its jacobians
|
|
J_n = J_(def_n)
|
|
J = J_(def_t)
|
|
# Define operators wrt to the reference domain
|
|
grad = lambda u: grad_(u, def_t)
|
|
sym_grad = lambda u: sym_grad_(u, def_t)
|
|
#inv_F = inv_F_(def_t)
|
|
div = lambda u: div_(u, def_t)
|
|
# Define domain velocity as a first order finite difference
|
|
w = (def_t - def_n)/tau
|
|
|
|
# All formulations contain GCL + Temam's stabilization term
|
|
if parameters['gcl_stab'] == True:
|
|
if parameters['jac_at_n'] == False:
|
|
print("Using time-derivative term approximated with Jacobian at n+1...")
|
|
F = fn.Constant(rho_f/tau)*J*fn.inner(u - u_n, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f/tau)*(J - J_n)*fn.inner(u, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f)*div(u_n - w)*fn.inner(u, v)*dx_f
|
|
else:
|
|
print("Using time-derivative term approximated with Jacobian at n...")
|
|
F = fn.Constant(rho_f/tau)*J_n*fn.inner(u - u_n, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f/tau)*(J - J_n)*fn.inner(u, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f)*div(u_n - w)*fn.inner(u, v)*dx_f
|
|
else:
|
|
raise Exception("Solver implemented only for cases with 'gcl_stab': True")
|
|
|
|
# Adding remaining term of the weak formulation
|
|
# The convective term is linearized!
|
|
F += fn.Constant(rho_f)*J*fn.dot(fn.dot(grad(u), u_n - w), v)*dx_f\
|
|
+ fn.Constant(2*mu)*J*fn.inner(sym_grad(u), sym_grad(v))*dx_f\
|
|
- p*div(v)*dx_f + q*div(u)*dx_f
|
|
|
|
# Add stabilization when required
|
|
if stabilization['type'] == None:
|
|
pass
|
|
else:
|
|
raise Exception("Other stabilization not required for this monolithic test model!")
|
|
|
|
# Obtain lhs and rhs
|
|
a, L = fn.lhs(F), fn.rhs(F)
|
|
# Assemble lhs of the formulation
|
|
A = fn.PETScMatrix()
|
|
fn.assemble(a, tensor = A)
|
|
# Define Krylov solver for the problem
|
|
if solver_to_use == 'LU':
|
|
print("Using LU Solver...")
|
|
solver = fn.LUSolver()
|
|
solver.set_operator(A)
|
|
elif solver_to_use == 'Krylov':
|
|
krylov_method = 'gmres'#'minres', 'gmres', 'bicgstab'
|
|
precond = 'icc'#'ilu', 'icc', 'hypre_amg', 'amg'
|
|
print("Using Krylov Solver --> Method: {}, Precond: {}\n".format(krylov_method, precond))
|
|
solver = fn.PETScKrylovSolver(krylov_method, precond)
|
|
solver.set_operator(A)
|
|
|
|
# TODO: Remove the zero-mean space for pressure fields
|
|
# Create normalized vector that spans the pressure null space (P)
|
|
p_ = fn.Function(P_)
|
|
null_p_vec = fn.Vector(p_.vector())
|
|
P_.dofmap().set(null_p_vec, 1.0)
|
|
null_p_vec *= 1.0/null_p_vec.norm("l2")
|
|
p_.vector()[:] = null_p_vec
|
|
# Create zero-vector that spans M by Function assigning it
|
|
null_func = fn.Function(M)
|
|
u_ = fn.Function(U_)
|
|
assigner_to_M.assign(null_func, [u_, p_])
|
|
null_vec = fn.Vector(null_func.vector())
|
|
# Create null space basis object and attach to PETSc matrix
|
|
null_space = fn.VectorSpaceBasis([null_vec])
|
|
fn.as_backend_type(A).set_nullspace(null_space)
|
|
|
|
fn.PETScOptions.set("ksp_monitor")
|
|
#solver.parameters['monitor_convergence'] = True
|
|
solver.parameters['maximum_iterations'] = 300
|
|
solver.parameters['relative_tolerance'] = 1E-2
|
|
solver.parameters['absolute_tolerance'] = 1E-4
|
|
|
|
else:
|
|
raise Exception("Solver type not implemented, available choices: 'LU', 'Krylov'")
|
|
|
|
# Define nicely the output file to save!
|
|
base_name = "iNSE_ALE"
|
|
base_name += "_{}_{}".format(def_name[0], def_name[1])
|
|
base_name += "_P{}P{}_GCL_{}".format(parameters['deg_U'], parameters['deg_P'], parameters['gcl_stab'])
|
|
|
|
if stabilization['type'] == None:
|
|
filename = base_name + "_NoStab.xdmf"
|
|
else:
|
|
raise Exception("Other stabilizations not implemented!")
|
|
# Create XDMF file (It saves multiple fields)
|
|
file_output = fn.XDMFFile('tests/rectangle/'+filename)
|
|
file_output.parameters['functions_share_mesh'] = True
|
|
|
|
# Starts from dt, we take at 0 the initial solution
|
|
print("Starting iteration over time...")
|
|
saved_defs = []
|
|
saved_vels = []
|
|
saved_pres = []
|
|
# Depending on time step tau, output some results to user!
|
|
if tau == 0.1 and T <= 1: nb = 2
|
|
else: nb = 1
|
|
# Saving velocity and deformation fields at time t = 0.
|
|
u_n, _ = aux_n.split(deepcopy=True)
|
|
file_output.write(u_n, 0.0)
|
|
energy_assembled = fn.assemble(J_n*fn.inner(u_n, u_n)*dx_f)*0.5*rho_f/tau
|
|
print("E (Energy) of initial velocity field: {:.5f} at time {:2.3f}\n".format(energy_assembled, 0.0))
|
|
aux_ = fn.Function(M.sub(0).collapse(), name="Deformation")
|
|
aux_.assign(fn.project(def_n, M.sub(0).collapse()))
|
|
file_output.write(aux_, 0.0)
|
|
# Main iteration
|
|
for count, t in enumerate(np.arange(tau, T+tau, tau)):
|
|
# re-Assemble lhs vector
|
|
fn.assemble(a, tensor=A)
|
|
[bc.apply(A) for bc in bdries]
|
|
# Assemble rhs vector (changes over time)
|
|
b = fn.assemble(L)
|
|
|
|
# if Krylov solver is used, then ortogonalize also b
|
|
if solver_to_use == 'Krylov':
|
|
null_space.orthogonalize(b)
|
|
|
|
[bc.apply(b) for bc in bdries]
|
|
|
|
# Solve linear system and save solution
|
|
aux = fn.Function(M, name="MixedSolution")
|
|
num_iters = solver.solve(aux.vector(), b)
|
|
# Split and save velocity field
|
|
u, p = aux.split(deepcopy=True)
|
|
file_output.write(u, t)
|
|
# Saving computed
|
|
aux_ = fn.Function(M.sub(0).collapse(), name="Deformation")
|
|
aux_.assign(fn.project(def_t, M.sub(0).collapse()))
|
|
file_output.write(aux_, t)
|
|
# Print only a couple of solutions to the user!
|
|
if count % nb == 0:
|
|
J_assembled = fn.assemble(J*dx_f)
|
|
J_n_assembled = fn.assemble(J_n*dx_f)
|
|
def_assembled = fn.assemble(fn.dot(def_t, def_t)*dx_f)
|
|
#w_assembled = fn.assemble(fn.dot(w, w)*dx_f)
|
|
energy_assembled = fn.assemble(J*fn.inner(u, u)*dx_f)*0.5*rho_f/tau
|
|
energy_n_assembled = fn.assemble(J_n*fn.inner(u_n, u_n)*dx_f)*0.5*rho_f/tau
|
|
energy_strain_assembled = fn.assemble(J*fn.inner(sym_grad(u), sym_grad(u))*dx_f)*2*mu
|
|
if parameters['jac_at_n'] == False:
|
|
energy_var_assembled = fn.assemble(J*fn.inner(u - u_n, u - u_n)*dx_f)*0.5*rho_f/tau
|
|
else:
|
|
energy_var_assembled = fn.assemble(J_n*fn.inner(u - u_n, u - u_n)*dx_f)*0.5*rho_f/tau
|
|
delta_value = energy_assembled - energy_n_assembled + energy_strain_assembled + energy_var_assembled
|
|
# Give info to user
|
|
print("J integrated over the domain: {:.3f} at time: {:2.3f}".format(J_assembled, t))
|
|
print("J_n integrated over the domain: {:.3f} at time: {:2.3f}".format(J_n_assembled, t))
|
|
print("Well-posedness condition 3J - J_n > 0: {} at time {:2.3f}".format(3*J_assembled - J_n_assembled > 0, t))
|
|
print("d (def. map) integrated over the domain: {:.3f} at time: {:2.3f}".format(def_assembled, t))
|
|
#print("w (vel. map)integrated over the domain: {:.3f} at time: {:2.3f}".format(w_assembled, t))
|
|
print("Strain energy E( eps(|u|^2) ): {:.5f} at time {:2.3f}".format(energy_strain_assembled, t))
|
|
print("Residual of scheme -> 'delta' value: {:.5f} at time {:2.3f}\n".format(delta_value, t))
|
|
# Saving results for postprocessing
|
|
if solver_to_use == 'LU':
|
|
saved_vels.append((u, t, delta_value , energy_strain_assembled, energy_var_assembled))
|
|
else:
|
|
saved_vels.append((u, t, delta_value , energy_strain_assembled, energy_var_assembled, num_iters))
|
|
|
|
# Update velocity profile
|
|
aux_n.assign(aux)
|
|
u_n, _ = fn.split(aux_n)
|
|
# Updates deformation mapping
|
|
def_n.t = t
|
|
def_t.t = t+tau
|
|
|
|
# Number of elements used and total time
|
|
print("Number of degrees of freedom: {}".format(M.dim()))
|
|
print("Solution Computed in {:.5} seconds".format(time.time()-start))
|
|
return saved_defs, saved_vels, saved_pres
|
|
|
|
|
|
def iNSE_CT_Solver(parameters, deformation, stabilization, solver_to_use='LU'):
|
|
""" This solver computes the incompressible Navier-Stokes (iNSE) problem in the ALE formalism.
|
|
To do so, the FEniCS framework is used and the convective term within iNSE is linealized.
|
|
To solve the problem, a Chorin-Temam scheme is employed with finite diferences for time
|
|
and FE for space.
|
|
The Chorin-Temam Scheme here is taken in the order: PPS_0 --> FVS_0 --> PPS_1 --> FVS_1...
|
|
...
|
|
|
|
Input Parameters
|
|
----------------
|
|
parameter: dict
|
|
dictionary containing several options for the scheme used.
|
|
deformation: dict
|
|
dictionary containing deformation mapping properties and options.
|
|
stabilization: dict
|
|
dictionary specifying stabilization used.
|
|
|
|
|
|
Output
|
|
----------------
|
|
None.
|
|
The solver print the results to console, and returns None, saving solutions and figures.
|
|
"""
|
|
# Starting time
|
|
start = time.time()
|
|
# Default parameters
|
|
mu, rho_f = parameters['mu'], parameters['rho_f']
|
|
Nx, Ny = parameters['Nx'], parameters['Ny']
|
|
tau, T = parameters['tau'], parameters['T']
|
|
|
|
def def_mesh_tuple(def_params: dict) -> tuple:
|
|
"""
|
|
Computes a list (len 2) containing the parametrized mesh velocity.
|
|
...
|
|
|
|
Input Parameters
|
|
----------------
|
|
def_params: dict
|
|
dictionary with keys 'expression' and 'output_name'.
|
|
|
|
"""
|
|
if 'expression' in def_params.keys():
|
|
if 'output_name' in def_params.keys():
|
|
output_name = def_params['output_name']
|
|
else:
|
|
output_name = ('X', 'Y')
|
|
return def_params['expression'], output_name
|
|
else:
|
|
raise Exception("Key 'expression' was not found in the dictionary!")
|
|
|
|
def exps_at_bdry(bdry_params: dict, def_params: dict) -> tuple:
|
|
"""
|
|
Computes the boundary FEniCS expressions for this particular Dirichlet Homogeneous case.
|
|
"""
|
|
u_inlet = fn.Constant((0.0, 0.0))
|
|
u_outlet = fn.Constant((0.0, 0.0))
|
|
u_bottom = fn.Constant((0.0, 0.0))
|
|
u_top = fn.Constant((0.0, 0.0))
|
|
return u_inlet, u_outlet, u_bottom, u_top
|
|
|
|
def corrected_vel(u, ut, p) -> fn.Function:
|
|
"""
|
|
Computes the corrected velocity field within CT update scheme.
|
|
"""
|
|
grad_p = fn.project(grad(p), U)
|
|
u.vector()[:] = ut.vector()[:] - (tau/rho_f)*grad_p.vector()[:]
|
|
return u
|
|
|
|
# Define boundary domain classes
|
|
class inlet(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[0], 0) and on_boundary
|
|
class outlet(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[0], 6) and on_boundary
|
|
class top(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[1], 1) and on_boundary
|
|
class bottom(fn.SubDomain):
|
|
def inside(self, x, on_boundary):
|
|
return fn.near(x[1], -1) and on_boundary
|
|
|
|
# Define mesh object
|
|
mesh_f = fn.RectangleMesh(fn.Point(0, -1), fn.Point(6, 1), Nx, Ny)
|
|
#h = fn.CellDiameter(mesh_f)
|
|
|
|
# Obtain mesh deformation tuple and output name
|
|
def_tuple, def_name = def_mesh_tuple(deformation)
|
|
|
|
# Define displacement expressions at previous and actual time steps
|
|
def_n = fn.Expression(def_tuple, t=0.0, domain=mesh_f, degree=2)
|
|
def_t = fn.Expression(def_tuple, t=tau, domain=mesh_f, degree=2)
|
|
|
|
# Functional Spaces
|
|
U = fn.VectorFunctionSpace(mesh_f, "CG", parameters['deg_U'])
|
|
P = fn.FunctionSpace(mesh_f, "CG", parameters['deg_P'])
|
|
|
|
# Boundary expressions in the reference domain
|
|
u_inflow, u_outflow, u_bottom, u_top = exps_at_bdry(parameters, deformation)
|
|
# Define and mark boundaries
|
|
bdry_markers = fn.MeshFunction('size_t', mesh_f, mesh_f.topology().dim()-1)
|
|
inlet().mark(bdry_markers, 1) # Inlet marked with label 1
|
|
outlet().mark(bdry_markers, 2) # Outlet marked with label 2
|
|
bottom().mark(bdry_markers, 3) # Bottom marked with label 3
|
|
top().mark(bdry_markers, 4) # Top marked with label 4
|
|
|
|
# Define measures
|
|
dx_f = fn.Measure('dx', domain=mesh_f, subdomain_data=bdry_markers)
|
|
#ds_f = fn.Measure('ds', domain=mesh_f, subdomain_data=bdry_markers)
|
|
|
|
# Define boundary conditions with pressure condition to obtain uniqueness (case LU solver!)
|
|
bdries_u = [
|
|
fn.DirichletBC(U, u_inflow, bdry_markers, 1),
|
|
fn.DirichletBC(U, u_outflow, bdry_markers, 2),
|
|
fn.DirichletBC(U, u_bottom, bdry_markers, 3),
|
|
fn.DirichletBC(U, u_top, bdry_markers, 4)
|
|
]
|
|
bdries_p = [
|
|
fn.DirichletBC(P, fn.Constant(0.), "near(x[0], 0) && near(x[1], 0)", 'pointwise')
|
|
]
|
|
|
|
# Trial/Test function definition
|
|
ut, v = fn.TrialFunction(U), fn.TestFunction(U)
|
|
p, q = fn.TrialFunction(P), fn.TestFunction(P)
|
|
ut_n, p_n = fn.Function(U, name='T. Velocity'), fn.Function(P, name='Pressure')
|
|
#u, u_n = fn.Function(U, name='C. Velocity'), fn.Function(U, name='C. Velocity')
|
|
# Define initial condition (non-zero)
|
|
if parameters["init_cond"] == True:
|
|
print("Using a non-zero initial condition (2D-Gaussian) for the velocity field!")
|
|
cpp_aux = '0.001*(1 - pow(x[1], 2))*x[0]*(6 - x[0])'
|
|
ut_n.assign(fn.Expression((cpp_aux, '0.'), domain=mesh_f, degree=2))
|
|
else:
|
|
pass
|
|
|
|
# Define normal for backflow stabilization term
|
|
#N = fn.FacetNormal(mesh_f)
|
|
# Compute its jacobians
|
|
J_n = J_(def_n)
|
|
J = J_(def_t)
|
|
# Define operators wrt to the reference domain
|
|
grad = lambda u: grad_(u, def_t)
|
|
grad_n = lambda u: grad_(u, def_n)
|
|
sym_grad = lambda u: sym_grad_(u, def_t)
|
|
div = lambda u: div_(u, def_t)
|
|
div_n = lambda u: div_(u, def_n)
|
|
# Define domain velocity as a first order finite difference
|
|
w = (def_t - def_n)/tau
|
|
|
|
# (FVS) Step. Finding tentative velocity field
|
|
# The formulation has GCL+Temam Stabilization term
|
|
print("Defining forms associated to Chorin-Temam scheme...")
|
|
if parameters['gcl_stab'] == True:
|
|
if parameters['jac_at_n'] == False:
|
|
F_fvs = fn.Constant(rho_f/tau)*J*fn.inner(ut - ut_n, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f/tau)*(J - J_n)*fn.inner(ut, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f)*div(ut_n - w)*fn.inner(ut, v)*dx_f
|
|
else:
|
|
F_fvs = fn.Constant(rho_f/tau)*J_n*fn.inner(ut - ut_n, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f/tau)*(J - J_n)*fn.inner(ut, v)*dx_f\
|
|
+ fn.Constant(0.5*rho_f)*div(ut_n - w)*fn.inner(ut, v)*dx_f
|
|
else:
|
|
raise Exception("Solver implemented only for cases with 'gcl_stab': True")
|
|
|
|
# Adding remaining term of the weak formulation
|
|
# The convective term is linearized!
|
|
F_fvs += fn.Constant(rho_f)*J*fn.dot(fn.dot(grad(ut), ut_n - w), v)*dx_f\
|
|
+ fn.Constant(2*mu)*J*fn.inner(sym_grad(ut), sym_grad(v))*dx_f\
|
|
- p_n*div_n(v)*dx_f
|
|
|
|
# Add stabilization when required
|
|
if stabilization['type'] == None:
|
|
pass
|
|
else:
|
|
raise Exception("Other stabilization not required for this CT test scheme!")
|
|
|
|
# (PPS) Step. Projecting into the div-free solution space
|
|
ut_pps = fn.Function(U) # Define tentative vel. object for pps form
|
|
F_pps = fn.Constant(tau/rho_f)*J_n*fn.inner(grad_n(p), grad_n(q))*dx_f\
|
|
+ div_n(ut_pps)*q*dx_f
|
|
|
|
# Obtain lhs and rhs
|
|
a_fvs, L_fvs = fn.lhs(F_fvs), fn.rhs(F_fvs)
|
|
a_pps, L_pps = fn.lhs(F_pps), fn.rhs(F_pps)
|
|
# Assemble lhs of the formulation
|
|
A_fvs = fn.PETScMatrix()
|
|
A_pps = fn.PETScMatrix()
|
|
fn.assemble(a_fvs, tensor = A_fvs)
|
|
fn.assemble(a_pps, tensor = A_pps)
|
|
|
|
# Define Krylov solver for the problem
|
|
# With 'Krylov' using 'gmres'+'ilu' the solver works smoothly
|
|
if solver_to_use == 'LU':
|
|
print("Using LU Solver...")
|
|
solver_fvs = fn.LUSolver()
|
|
solver_pps = fn.LUSolver()
|
|
[solver.set_operator(A) for solver, A in zip([solver_fvs, solver_pps], [A_fvs, A_pps])]
|
|
elif solver_to_use == 'Krylov':
|
|
krylov_method = 'gmres'#'minres', 'gmres', 'bicgstab'
|
|
precond = 'amg'#'ilu', 'icc', 'hypre_amg', 'petsc_amg'
|
|
print("Using Krylov Solver --> Method: {}, Precond: {}\n".format(krylov_method, precond))
|
|
solver_fvs = fn.PETScKrylovSolver(krylov_method, precond)
|
|
solver_pps = fn.PETScKrylovSolver(krylov_method, precond)
|
|
for solver, A in zip([solver_fvs, solver_pps], [A_fvs, A_pps]):
|
|
solver.set_operator(A)
|
|
solver.parameters['report'] = True
|
|
#solver.parameters['monitor_convergence'] = True
|
|
solver.parameters['maximum_iterations'] = 2000
|
|
solver.parameters['relative_tolerance'] = 1E-4
|
|
solver.parameters['absolute_tolerance'] = 1E-7
|
|
#solver.parameters.str(True)
|
|
else:
|
|
raise Exception("Solver type not implemented, available choices: 'LU', 'Krylov'")
|
|
|
|
# Define nicely the output file to save!
|
|
base_name = "iNSE_ALE_CT"
|
|
base_name += "_{}_{}".format(def_name[0], def_name[1])
|
|
base_name += "_P{}P{}_GCL_{}".format(parameters['deg_U'], parameters['deg_P'], parameters['gcl_stab'])
|
|
|
|
if stabilization['type'] == None:
|
|
filename = base_name + "_NoStab.xdmf"
|
|
else:
|
|
raise Exception("Other stabilizations not implemented!")
|
|
# Create XDMF file (It saves multiple fields)
|
|
file_output = fn.XDMFFile('tests/rectangle/'+filename)
|
|
file_output.parameters['functions_share_mesh'] = True
|
|
|
|
# Starts from dt, we take at 0 the initial solution
|
|
print("Starting iteration over time...")
|
|
saved_defs = []
|
|
saved_vels = []
|
|
saved_pres = []
|
|
# Depending on time step tau, save and export some solutions
|
|
if tau == 0.1 and T <= 1: nb = 2
|
|
else: nb = 1
|
|
|
|
##-----Saving deformation field at time t = 0-----##
|
|
aux_ = fn.Function(U, name="Deformation")
|
|
aux_.assign(fn.project(def_n, U))
|
|
file_output.write(aux_, 0.0)
|
|
# Saving velocity (t., div-free) and pressure fields
|
|
# by solving an initial PPS
|
|
fn.assemble(a_pps, tensor=A_pps)
|
|
[bc.apply(A_pps) for bc in bdries_p]
|
|
# Assemble rhs (changes over time)
|
|
b_pps = fn.assemble(L_pps)
|
|
[bc.apply(b_pps) for bc in bdries_p]
|
|
# Solve linear system
|
|
p = fn.Function(P, name='Pressure')
|
|
num_iters_pps = solver_pps.solve(p.vector(), b_pps)
|
|
file_output.write(p, 0.0)
|
|
# Compute and save updated (t.) velocity
|
|
corrected_vel(ut_n, ut_n, p)
|
|
file_output.write(ut_n, 0.0)
|
|
#corrected_vel(u_n, ut_n, p_n)
|
|
#file_output.write(u_n, 0.0)
|
|
energy_assembled = fn.assemble(J_n*fn.inner(ut_n, ut_n)*dx_f)*0.5*rho_f/tau
|
|
print("E (Energy) of initial t. velocity: {:.5f} at time {:2.3f}\n".format(energy_assembled, 0.0))
|
|
## -------------------------------------------- ##
|
|
## ----------------Main iteration---------------##
|
|
for count, t in enumerate(np.arange(tau, T+tau, tau)):
|
|
# (PPS) Pressure-Projection Step
|
|
# Update t. velocity
|
|
ut_pps.assign(ut_n)
|
|
# Assemble lhs
|
|
fn.assemble(a_pps, tensor=A_pps)
|
|
[bc.apply(A_pps) for bc in bdries_p]
|
|
# Assemble rhs (changes over time)
|
|
b_pps = fn.assemble(L_pps)
|
|
[bc.apply(b_pps) for bc in bdries_p]
|
|
# Solve linear system
|
|
p = fn.Function(P, name='Pressure')
|
|
num_iters_pps = solver_pps.solve(p.vector(), b_pps)
|
|
# Save pressure
|
|
file_output.write(p, t)
|
|
|
|
# (FVS) Fluid-Viscous Step
|
|
# Update pressure
|
|
p_n.assign(p)
|
|
# re-Assemble lhs vector
|
|
fn.assemble(a_fvs, tensor=A_fvs)
|
|
[bc.apply(A_fvs) for bc in bdries_u]
|
|
# Assemble rhs vector (changes over time)
|
|
b_fvs = fn.assemble(L_fvs)
|
|
[bc.apply(b_fvs) for bc in bdries_u]
|
|
# Solve linear system and save solution
|
|
ut = fn.Function(U, name='T. Velocity')
|
|
num_iters_fvs = solver_fvs.solve(ut.vector(), b_fvs)
|
|
|
|
# (CVS) Computation Corrected Velocity (Optional)
|
|
#corrected_vel(u, ut, p)
|
|
|
|
# Save deformation (proj. into U)
|
|
aux_ = fn.Function(U, name="Deformation")
|
|
aux_.assign(fn.project(def_t, U))
|
|
file_output.write(aux_, t)
|
|
# Save t. velocity
|
|
file_output.write(ut, t)
|
|
#file_output.write(u, t)
|
|
|
|
# Print and save solutions to the user!
|
|
if count % nb == 0:
|
|
J_assembled = fn.assemble(J*dx_f)
|
|
J_n_assembled = fn.assemble(J_n*dx_f)
|
|
def_assembled = fn.assemble(fn.dot(def_t, def_t)*dx_f)
|
|
#w_assembled = fn.assemble(fn.dot(w, w)*dx_f)
|
|
energy_assembled = fn.assemble(J*fn.inner(ut, ut)*dx_f)*0.5*rho_f/tau
|
|
energy_n_assembled = fn.assemble(J_n*fn.inner(ut_n, ut_n)*dx_f)*0.5*rho_f/tau
|
|
energy_strain_assembled = fn.assemble(J*fn.inner(sym_grad(ut), sym_grad(ut))*dx_f)*2*mu
|
|
pres_assembled = fn.assemble(J_n*fn.inner(grad_n(p), grad_n(p))*dx_f)*0.5*tau/rho_f
|
|
delta_value = energy_assembled - energy_n_assembled + energy_strain_assembled + pres_assembled
|
|
if parameters['jac_at_n'] == False:
|
|
energy_var_assembled = fn.assemble(J*fn.inner(ut - ut_n, ut - ut_n)*dx_f)*0.5*rho_f/tau
|
|
else:
|
|
energy_var_assembled = fn.assemble(J_n*fn.inner(ut - ut_n, ut - ut_n)*dx_f)*0.5*rho_f/tau
|
|
# Providing info to user
|
|
print("J integrated over the domain: {:.3f} at time: {:2.3f}".format(J_assembled, t))
|
|
print("J_n integrated over the domain: {:.3f} at time: {:2.3f}".format(J_n_assembled, t))
|
|
print("Well-posedness condition 3J - J_n > 0: {} at time {:2.3f}".format(3*J_assembled - J_n_assembled > 0, t))
|
|
print("d (def. map) integrated over the domain: {:.3f} at time: {:2.3f}".format(def_assembled, t))
|
|
#print("w (vel. map)integrated over the domain: {:.3f} at time: {:2.3f}".format(w_assembled, t))
|
|
print("Strain energy E( eps(|u|^2) ): {:.5f} at time {:2.3f}".format(energy_strain_assembled, t))
|
|
print("Residual of scheme -> 'delta' value: {:.5f} at time {:2.3f}\n".format(delta_value, t))
|
|
# Saving results
|
|
if solver_to_use == 'LU':
|
|
saved_vels.append((ut, t, delta_value, energy_strain_assembled, energy_var_assembled))
|
|
else:
|
|
saved_vels.append((ut, t, delta_value, energy_strain_assembled, energy_var_assembled, num_iters_fvs, num_iters_pps))
|
|
|
|
# Update t. velocity
|
|
ut_n.assign(ut)
|
|
#u_n.assign(u)
|
|
# Updates deformation mapping
|
|
def_n.t = t
|
|
def_t.t = t+tau
|
|
|
|
# Number of elements used and total time
|
|
print("DoF for Velocity Space (U): {}".format(U.dim()))
|
|
print("DoF for Pressure Space (P): {}".format(P.dim()))
|
|
print("Solution Computed in {:.5} seconds".format(time.time()-start))
|
|
return saved_defs, saved_vels, saved_pres
|
|
|
|
def testing_comparative_plots(dict_fields_mt: dict, dict_fields_ct: dict, *args) -> None:
|
|
"""
|
|
Computes error figures between the schemes with Jacobian evaluation at times {n, n+1}, i.e.
|
|
with explicit and implicit treatment of the time-derivative term.
|
|
...
|
|
|
|
Input Parameters
|
|
----------------
|
|
dict_fields_mt : tuple -> (dict, dict, dict)
|
|
tuple containing dictionaries of fields (def, vel, pre) from the iNSE ALE Monolithic solver.
|
|
|
|
dict_fields_ct : tupel -> (dict, dict, dict)
|
|
tuple containing dictionaries of fields (def, vel, pre) from the iNSE ALE CT solver.
|
|
"""
|
|
print("\nComputing error figures for both schemes...\n")
|
|
# Creates dictionaries for "expanding" and "contracting" functions
|
|
# Time dictionary
|
|
time_list = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
# Numerical dissipation dictionary --> int (| J^{n+1} u^{n+1} - J^{n} u^{n} |)
|
|
dis_value_exp = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
dis_value_con = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
# Physical dissipation dictionary --> int ( 2*mu* | eps(u) |^2 )
|
|
eps_value_exp = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
eps_value_con = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
# 'delta' values dictionary
|
|
delta_value_exp = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
delta_value_con = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
# Normalized 'delta' values dictionary (w.r.t physical dissipation)
|
|
delta_hat_value_exp = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
delta_hat_value_con = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
# Normalized 'delta' values dictionary (w.r.t numerical dissipation)
|
|
delta_hat_dis_exp = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
delta_hat_dis_con = {
|
|
'MT': {jac_name: [] for jac_name in dict_fields_mt.keys()},
|
|
'CT': {jac_name: [] for jac_name in dict_fields_ct.keys()}
|
|
}
|
|
# Iterate over the different cases in MT (Monolithic) and CT (Chorin-Temam) schemes
|
|
for scheme in ['MT', 'CT']: # 'MT', 'CT'
|
|
if scheme == 'MT':
|
|
dict_fields = dict_fields_mt
|
|
else:
|
|
dict_fields = dict_fields_ct
|
|
# Iterate over different cases 'j_implicit', 'j_explicit'
|
|
for jac_name in dict_fields.keys():
|
|
dict_vel = dict_fields[jac_name][1]
|
|
for i in range(len(dict_vel['High_Osc_exp'])):
|
|
# Saving time steps
|
|
time_list[scheme][jac_name].append(dict_vel['High_Osc_exp'][i][1])
|
|
# Add 'delta' values
|
|
delta_value_exp[scheme][jac_name].append(dict_vel['High_Osc_exp'][i][2])
|
|
delta_value_con[scheme][jac_name].append(dict_vel['High_Osc_con'][i][2])
|
|
# Add physical dissipation
|
|
eps_value_exp[scheme][jac_name].append(dict_vel['High_Osc_exp'][i][3])
|
|
eps_value_con[scheme][jac_name].append(dict_vel['High_Osc_con'][i][3])
|
|
# Add numerical dissipation
|
|
dis_value_exp[scheme][jac_name].append(dict_vel['High_Osc_exp'][i][4])
|
|
dis_value_con[scheme][jac_name].append(dict_vel['High_Osc_con'][i][4])
|
|
# Compute normalized 'delta' values w.r.t the physical dissipation
|
|
delta_hat_value_exp[scheme][jac_name].append(dict_vel['High_Osc_exp'][i][2]/dict_vel['High_Osc_exp'][i][3])
|
|
delta_hat_value_con[scheme][jac_name].append(dict_vel['High_Osc_con'][i][2]/dict_vel['High_Osc_con'][i][3])
|
|
# Compute normalized 'delta' values w.r.t the numerical dissipation
|
|
delta_hat_dis_exp[scheme][jac_name].append(dict_vel['High_Osc_exp'][i][2]/dict_vel['High_Osc_exp'][i][4])
|
|
delta_hat_dis_con[scheme][jac_name].append(dict_vel['High_Osc_con'][i][2]/dict_vel['High_Osc_con'][i][4])
|
|
|
|
# Auxiliary indexes for plotting
|
|
# TODO: Test this line and automatize the "100"!
|
|
init = 0
|
|
out = 101 #len(time_list['CT']['j_implicit']) # at tau = 0.01, T = 2 --> 200 time steps
|
|
|
|
# ---------------------- #
|
|
# DELTA VALUE CASE
|
|
# ---------------------- #
|
|
# CASE CT
|
|
plt.subplots(figsize=(10, 4))
|
|
# Top
|
|
lf_ax = plt.subplot(1, 2, 2)
|
|
lf_j_exp, = lf_ax.plot(time_list['CT']['j_explicit'][init:out], delta_value_exp['CT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['CT']['j_implicit'][init:out], delta_value_exp['CT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"$\delta_{CT}$ def. I", color='blue')
|
|
lf_ax.set_title(r"Residual Energy $\delta_{CT}$")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['CT']['j_explicit'][init:out], delta_value_con['CT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['CT']['j_implicit'][init:out], delta_value_con['CT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"$\delta_{CT}$ def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'CT $\star\star = n$', r'CT $\star\star = n+1$', r'CT $\star\star = n$', r'CT $\star\star = n+1$'), loc='center right')
|
|
plt.grid(True)
|
|
# CASE M
|
|
# Bottom
|
|
lf_ax = plt.subplot(1, 2, 1)
|
|
lf_j_exp, = lf_ax.plot(time_list['MT']['j_explicit'][init:out], delta_value_exp['MT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['MT']['j_implicit'][init:out], delta_value_exp['MT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"$\delta_{M}$ def. I", color='blue')
|
|
lf_ax.set_title(r"Residual Energy $\delta_{M}$")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['MT']['j_explicit'][init:out], delta_value_con['MT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['MT']['j_implicit'][init:out], delta_value_con['MT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"$\delta_{M}$ def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'M $\star\star = n$', r'M $\star\star = n+1$', r'M $\star\star = n$', r'M $\star\star = n+1$'), loc='center right')
|
|
plt.grid(True)
|
|
plt.subplots_adjust(hspace=0.8)
|
|
plt.tight_layout()
|
|
plt.savefig('Comparison_Delta_Value_GCL_{}_Solver_{}.png'.format(args[0], args[1]), dpi=300)
|
|
plt.close()
|
|
|
|
# ---------------------- #
|
|
# NORMALIZED DELTA CASE (w.r.t physical dissipation)
|
|
# ---------------------- #
|
|
plt.subplots(figsize=(10, 4))
|
|
# CASE CT
|
|
# Top
|
|
lf_ax = plt.subplot(1, 2, 2)
|
|
lf_j_exp, = lf_ax.plot(time_list['CT']['j_explicit'][init:out], delta_hat_value_exp['CT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['CT']['j_implicit'][init:out], delta_hat_value_exp['CT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"$\hat{\delta}_{CT}$ def. I", color='blue')
|
|
lf_ax.set_title(r"Normalized value $\hat{\delta}_{CT}$")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['CT']['j_explicit'][init:out], delta_hat_value_con['CT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['CT']['j_implicit'][init:out], delta_hat_value_con['CT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"$\hat{\delta}_{CT}$ def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'CT $\star\star = n$', r'CT $\star\star = n+1$', r'CT $\star\star = n$', r'CT $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
# CASE M
|
|
# Bottom
|
|
lf_ax = plt.subplot(1, 2, 1)
|
|
lf_j_exp, = lf_ax.plot(time_list['MT']['j_explicit'][init:out], delta_hat_value_exp['MT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['MT']['j_implicit'][init:out], delta_hat_value_exp['MT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"$\hat{\delta}_{M}$ def. I", color='blue')
|
|
lf_ax.set_title(r"Normalized value $\hat{\delta}_{M}$")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['MT']['j_explicit'][init:out], delta_hat_value_con['MT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['MT']['j_implicit'][init:out], delta_hat_value_con['MT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"$\hat{\delta}_{M}$ def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'M $\star\star = n$', r'M $\star\star = n+1$', r'M $\star\star = n$', r'M $\star\star = n+1$'), loc='upper right')
|
|
plt.grid(True)
|
|
plt.subplots_adjust(hspace=0.8)
|
|
plt.tight_layout()
|
|
plt.savefig('Comparison_Delta_Hat_Value_GCL_{}_solver_{}.png'.format(args[0], args[1]), dpi=300)
|
|
plt.close()
|
|
|
|
# ---------------------- #
|
|
# NORMALIZED DELTA CASE (w.r.t numerical dissipation)
|
|
# ---------------------- #
|
|
plt.subplots(figsize=(10, 4))
|
|
# CASE CT
|
|
# Top
|
|
lf_ax = plt.subplot(1, 2, 2)
|
|
lf_j_exp, = lf_ax.plot(time_list['CT']['j_explicit'][init:out], delta_hat_dis_exp['CT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['CT']['j_implicit'][init:out], delta_hat_dis_exp['CT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"$\hat{\delta}_{CT}$ dis. def. I", color='blue')
|
|
lf_ax.set_title(r"Normalized value $\hat{\delta}_{CT}$ (dis)")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['CT']['j_explicit'][init:out], delta_hat_dis_con['CT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['CT']['j_implicit'][init:out], delta_hat_dis_con['CT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"$\hat{\delta}_{CT}$ def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'CT $\star\star = n$', r'CT $\star\star = n+1$', r'CT $\star\star = n$', r'CT $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
# CASE M
|
|
# Bottom
|
|
lf_ax = plt.subplot(1, 2, 1)
|
|
lf_j_exp, = lf_ax.plot(time_list['MT']['j_explicit'][init:out], delta_hat_dis_exp['MT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['MT']['j_implicit'][init:out], delta_hat_dis_exp['MT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"$\hat{\delta}_{M}$ dis. def. I", color='blue')
|
|
lf_ax.set_title(r"Normalized value $\hat{\delta}_{M}$ (dis)")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['MT']['j_explicit'][init:out], delta_hat_dis_con['MT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['MT']['j_implicit'][init:out], delta_hat_dis_con['MT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"$\hat{\delta}_{M}$ def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'M $\star\star = n$', r'M $\star\star = n+1$', r'M $\star\star = n$', r'M $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
plt.subplots_adjust(hspace=0.8)
|
|
plt.tight_layout()
|
|
plt.savefig('Comparison_Delta_Hat_Dis_Value_GCL_{}_solver_{}.png'.format(args[0], args[1]), dpi=300)
|
|
plt.close()
|
|
|
|
# ---------------------- #
|
|
# PHYSICAL DISSIPATION PLOT
|
|
# ---------------------- #
|
|
plt.subplots(figsize=(10, 4))
|
|
# CASE CT
|
|
# Top
|
|
lf_ax = plt.subplot(1, 2, 2)
|
|
lf_j_exp, = lf_ax.plot(time_list['CT']['j_explicit'][init:out], eps_value_exp['CT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['CT']['j_implicit'][init:out], eps_value_exp['CT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"On def. I", color='blue')
|
|
lf_ax.set_title(r"Physical Dissipation (CT)")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['CT']['j_explicit'][init:out], eps_value_con['CT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['CT']['j_implicit'][init:out], eps_value_con['CT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"On def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'CT $\star\star = n$', r'CT $\star\star = n+1$', r'CT $\star\star = n$', r'CT $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
# CASE M
|
|
# Bottom
|
|
lf_ax = plt.subplot(1, 2, 1)
|
|
lf_j_exp, = lf_ax.plot(time_list['MT']['j_explicit'][init:out], eps_value_exp['MT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['MT']['j_implicit'][init:out], eps_value_exp['MT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"On dis. def. I", color='blue')
|
|
lf_ax.set_title(r"Physical Dissipation (M)")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['MT']['j_explicit'][init:out], eps_value_con['MT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['MT']['j_implicit'][init:out], eps_value_con['MT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"On def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'M $\star\star = n$', r'M $\star\star = n+1$', r'M $\star\star = n$', r'M $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
plt.subplots_adjust(hspace=0.8)
|
|
plt.tight_layout()
|
|
plt.savefig('Comparison_Eps_Value_GCL_{}_solver_{}.png'.format(args[0], args[1]), dpi=300)
|
|
plt.close()
|
|
|
|
# ---------------------- #
|
|
# NUMERICAL DISSIPATION PLOT
|
|
# ---------------------- #
|
|
plt.subplots(figsize=(10, 4))
|
|
# CASE CT
|
|
# Top
|
|
lf_ax = plt.subplot(1, 2, 2)
|
|
lf_j_exp, = lf_ax.plot(time_list['CT']['j_explicit'][init:out], dis_value_exp['CT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['CT']['j_implicit'][init:out], dis_value_exp['CT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"On def. I", color='blue')
|
|
lf_ax.set_title(r"Numerical Dissipation (CT)")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['CT']['j_explicit'][init:out], dis_value_con['CT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['CT']['j_implicit'][init:out], dis_value_con['CT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"On def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'CT $\star\star = n$', r'CT $\star\star = n+1$', r'CT $\star\star = n$', r'CT $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
# CASE M
|
|
# Bottom
|
|
lf_ax = plt.subplot(1, 2, 1)
|
|
lf_j_exp, = lf_ax.plot(time_list['MT']['j_explicit'][init:out], dis_value_exp['MT']['j_explicit'][init:out], color='blue', ls='--')
|
|
lf_j_imp, = lf_ax.plot(time_list['MT']['j_implicit'][init:out], dis_value_exp['MT']['j_implicit'][init:out], color='blue')
|
|
lf_ax.set_xlabel(r"Simulation time [t]")
|
|
lf_ax.set_ylabel(r"On def. I", color='blue')
|
|
lf_ax.set_title(r"Numerical Dissipation (M)")
|
|
lf_ax.tick_params(axis='y', labelcolor='blue')
|
|
lf_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
rh_ax = lf_ax.twinx()
|
|
rh_j_exp, = rh_ax.plot(time_list['MT']['j_explicit'][init:out], dis_value_con['MT']['j_explicit'][init:out], color='red', ls='--')
|
|
rh_j_imp, = rh_ax.plot(time_list['MT']['j_implicit'][init:out], dis_value_con['MT']['j_implicit'][init:out], color='red')
|
|
rh_ax.set_ylabel(r"On def. II", color='red')
|
|
rh_ax.tick_params(axis='y', labelcolor='red')
|
|
rh_ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
|
|
plt.legend((lf_j_exp, lf_j_imp, rh_j_exp, rh_j_imp), \
|
|
(r'M $\star\star = n$', r'M $\star\star = n+1$', r'M $\star\star = n$', r'M $\star\star = n+1$'), loc='lower right')
|
|
plt.grid(True)
|
|
plt.subplots_adjust(hspace=0.8)
|
|
plt.tight_layout()
|
|
plt.savefig('Comparison_Num_Dis_Value_GCL_{}_solver_{}.png'.format(args[0], args[1]), dpi=300)
|
|
plt.close()
|
|
print("\nAll computations done and figures saved...\n")
|
|
print("\n------------------------------------------\n")
|
|
return None |