iNSE-ALE-Article/solvers.py

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