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