iNSE-ALE-Article/ALE Testing - Article.ipynb

In [ ]:

from main import testing_ALE, testing_defs_ALE, testing_defs_ALE_CT
#testing_defs_ALE(solver_to_use='LU')
#testing_defs_ALE_CT(solver_to_use='LU')
dict_mt, dict_ct = testing_ALE(solver_to_use='LU') # 'LU', 'Krylov'

In [ ]:

import matplotlib.pyplot as plt
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}$", 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), (r'CT $\bigstar\bigstar = n$', r'CT $\bigstar\bigstar = 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}$", 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), (r'M $\bigstar\bigstar = n$', r'M $\bigstar\bigstar = 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}$", 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), (r'CT $\bigstar\bigstar = n$', r'CT $\bigstar\bigstar = 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}$", 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), (r'M $\bigstar\bigstar = n$', r'M $\bigstar\bigstar = 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. ", 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), (r'CT $\bigstar\bigstar = n$', r'CT $\bigstar\bigstar = 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.", 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), (r'M $\bigstar\bigstar = n$', r'M $\bigstar\bigstar = 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"CT", 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), \
        (r'CT $\star\star = n$', r'CT $\bigstar\bigstar = 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"M", 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), (r'M $\bigstar\bigstar = n$', r'M $\bigstar\bigstar = 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"CT", 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), (r'CT $\bigstar\bigstar = n$', r'CT $\bigstar\bigstar = 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"M", 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), (r'M $\bigstar\bigstar = n$', r'M $\bigstar\bigstar = 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

testing_comparative_plots(dict_mt, dict_ct, True, 'LU')

In [ ]: