iNSE-ALE-Article/solvers.py

1005 lines
48 KiB
Python
Raw Normal View History

"""
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