From 23db005da728e58330cadb75fb24300a1a822762 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Wed, 10 Jul 2024 18:12:44 +0100 Subject: [PATCH 1/3] A series of useful python scripts based on xarray for plotting data from moment_kinetics simulations. --- xarray_post_processing/plot_error_data.py | 119 ++++++++ .../plot_integration_error_data.py | 62 ++++ .../plot_many_collisions.py | 267 ++++++++++++++++++ xarray_post_processing/plot_mk_utils.py | 148 ++++++++++ xarray_post_processing/plot_sd.py | 176 ++++++++++++ xarray_post_processing/plot_wall.py | 150 ++++++++++ xarray_post_processing/xarray_mk_utils.py | 33 +++ 7 files changed, 955 insertions(+) create mode 100644 xarray_post_processing/plot_error_data.py create mode 100644 xarray_post_processing/plot_integration_error_data.py create mode 100644 xarray_post_processing/plot_many_collisions.py create mode 100644 xarray_post_processing/plot_mk_utils.py create mode 100644 xarray_post_processing/plot_sd.py create mode 100644 xarray_post_processing/plot_wall.py create mode 100644 xarray_post_processing/xarray_mk_utils.py diff --git a/xarray_post_processing/plot_error_data.py b/xarray_post_processing/plot_error_data.py new file mode 100644 index 000000000..1dfac40f9 --- /dev/null +++ b/xarray_post_processing/plot_error_data.py @@ -0,0 +1,119 @@ +# -*- coding: utf-8 -*- +import matplotlib +# Force matplotlib to not use any Xwindows backend. +matplotlib.use('Agg') # this line allows plots to be made without using a display environment variable +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np +import h5py +from plot_mk_utils import plot_1d_list_pdf, plot_1d_loglog_list_pdf +from plot_mk_utils import plot_1d_semilog_list_pdf + + +def get_fkpl_error_data(filename): + f = h5py.File(filename,'r') + print(f.keys()) + ncore = np.copy(f['ncore'][...]) + ngrid = np.copy(f['ngrid'][...]) + print("ngrid: ",ngrid) + nelement_list = np.copy(f['nelement_list'][:]) + print("nelement_list: ",nelement_list[:]) + max_C_err = np.copy(f['max_C_err'][:]) + max_G_err = np.copy(f['max_G_err'][:]) + max_H_err = np.copy(f['max_H_err'][:]) + max_dHdvpa_err = np.copy(f['max_dHdvpa_err'][:]) + max_dHdvperp_err = np.copy(f['max_dHdvperp_err'][:]) + max_d2Gdvperpdvpa_err = np.copy(f['max_d2Gdvperpdvpa_err'][:]) + max_d2Gdvpa2_err = np.copy(f['max_d2Gdvpa2_err'][:]) + max_d2Gdvperp2_err = np.copy(f['max_d2Gdvperp2_err'][:]) + L2_C_err = np.copy(f['L2_C_err'][:]) + L2_G_err = np.copy(f['L2_G_err'][:]) + L2_H_err = np.copy(f['L2_H_err'][:]) + L2_dHdvpa_err = np.copy(f['L2_dHdvpa_err'][:]) + L2_dHdvperp_err = np.copy(f['L2_dHdvperp_err'][:]) + L2_d2Gdvperpdvpa_err = np.copy(f['L2_d2Gdvperpdvpa_err'][:]) + L2_d2Gdvpa2_err = np.copy(f['L2_d2Gdvpa2_err'][:]) + L2_d2Gdvperp2_err = np.copy(f['L2_d2Gdvperp2_err'][:]) + expected_diff = np.copy(f['expected_diff'][:]) + expected_integral = np.copy(f['expected_integral'][:]) + calculate_times = np.copy(f['calculate_times'][:]) + init_times = np.copy(f['init_times'][:]) + expected_t_2 = np.copy(f['expected_t_2'][:]) + expected_t_3 = np.copy(f['expected_t_3'][:]) + n_err = np.copy(f['n_err'][:]) + u_err = np.copy(f['u_err'][:]) + p_err = np.copy(f['p_err'][:]) + print(p_err) + + nelement_string = "N_{\\rm EL}" #\\scriptscriptstyle + ngrid_string = "N_{\\rm GR}" + + + file = filename+".plots.pdf" + pdf = PdfPages(file) + + marker_list = ['r--o','b--s','g--.','m--x','c--v','k'] + C_list = [max_C_err,L2_C_err,n_err,u_err,p_err,expected_diff] + nelements = [nelement_list for item in C_list] + ylab_list = ["$\\epsilon_{\\infty}(C[F,F])$", "$\\epsilon_{L_2}(C[F,F])$", + "$|\\Delta n|$", + "$|\\Delta u_{||}|$", + "$|\\Delta p |$", + "$(1/"+nelement_string+")^{"+ngrid_string+"-1}$"] + plot_1d_loglog_list_pdf (nelements,C_list,marker_list,"$"+ nelement_string+"$", pdf, + title='',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = nelement_list, yticks = None, + markersize=10, legend_title="", use_legend=True,loc_opt='lower left', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.05, 0.05), legend_fontsize=15, ncol_opt=1) + + marker_list = ['r--o','b--s','m-o','k-s'] + time_list = [init_times, calculate_times, expected_t_3, expected_t_2] + nelements = [nelement_list for item in time_list] + ylab_list = ["time/init (ms)", "time/step (ms)", + "$"+nelement_string+"^3$", + "$"+nelement_string+"^2$"] + plot_1d_loglog_list_pdf (nelements,time_list,marker_list,"$"+ nelement_string+"$", pdf, + title='',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = nelement_list, yticks = None, + markersize=10, legend_title="", use_legend=True,loc_opt='upper left', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.05, 0.95), legend_fontsize=15, ncol_opt=1) + + + marker_list = ['r--o','b--^','g--.','m--x','c--v','k'] + Infnorm_list = [max_dHdvpa_err, + max_dHdvperp_err,max_d2Gdvperpdvpa_err,max_d2Gdvpa2_err,max_d2Gdvperp2_err, + #expected_diff, + expected_integral] + nelements = [nelement_list for item in Infnorm_list] + ylab_list = ["$\\epsilon_{\\infty}(d H / d v_{||})$","$\\epsilon_{\\infty}(d H / d v_{\\perp})$", + "$\\epsilon_{\\infty}(d^2 G / d v_{\\perp} d v_{||})$", + "$\\epsilon_{\\infty}(d^2 G / d v^2_{||})$", + "$\\epsilon_{\\infty}(d^2 G / d v^2_{\\perp})$", + #"$(1/"+nelement_string+")^{"+ngrid_string+"-1}$", + "$(1/"+nelement_string+")^{"+ngrid_string+"+1}$"] + plot_1d_loglog_list_pdf (nelements,Infnorm_list,marker_list,"$"+ nelement_string+"$", pdf, + title='',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = nelement_list, yticks = None, + markersize=10, legend_title="", use_legend=True,loc_opt='lower left', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.05, 0.05), legend_fontsize=15, ncol_opt=1) + + L2norm_list = [L2_dHdvpa_err, + L2_dHdvperp_err,L2_d2Gdvperpdvpa_err,L2_d2Gdvpa2_err,L2_d2Gdvperp2_err, + #expected_diff, + expected_integral] + nelements = [nelement_list for item in L2norm_list] + ylab_list = ["$\\epsilon_{L_2}(d H / d v_{||})$","$\\epsilon_{L_2}(d H / d v_{\\perp})$", + "$\\epsilon_{L_2}(d^2 G / d v_{\\perp} d v_{||})$", + "$\\epsilon_{L_2}(d^2 G / d v^2_{||})$", + "$\\epsilon_{L_2}(d^2 G / d v^2_{\\perp})$", + #"$(1/"+nelement_string+")^{"+ngrid_string+"-1}$", + "$(1/"+nelement_string+")^{"+ngrid_string+"+1}$"] + plot_1d_loglog_list_pdf (nelements,L2norm_list,marker_list,"$"+ nelement_string+"$", pdf, + title='',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = nelement_list, yticks = None, + markersize=10, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) + pdf.close() + print(file) + + f.close() + return None + +workdir = "" +filename = workdir + "moment_kinetics_collisions/fkpl_error_data_ngrid_9_ncore_2.h5" +get_fkpl_error_data(filename) diff --git a/xarray_post_processing/plot_integration_error_data.py b/xarray_post_processing/plot_integration_error_data.py new file mode 100644 index 000000000..0ddd04f42 --- /dev/null +++ b/xarray_post_processing/plot_integration_error_data.py @@ -0,0 +1,62 @@ +# -*- coding: utf-8 -*- +import matplotlib +# Force matplotlib to not use any Xwindows backend. +matplotlib.use('Agg') # this line allows plots to be made without using a display environment variable +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np +import h5py +from plot_mk_utils import plot_1d_list_pdf, plot_1d_loglog_list_pdf +from plot_mk_utils import plot_1d_semilog_list_pdf + + +def get_fkpl_integration_error_data(filename): + f = h5py.File(filename,'r') + print(f.keys()) + ncore = np.copy(f['ncore'][...]) + ngrid = np.copy(f['ngrid'][...]) + print("ngrid: ",ngrid) + nelement_list = np.copy(f['nelement_list'][:]) + print("nelement_list: ",nelement_list[:]) + max_dHdvpa_err = np.copy(f['max_dHdvpa_err'][:]) + max_dHdvperp_err = np.copy(f['max_dHdvperp_err'][:]) + max_d2Gdvperpdvpa_err = np.copy(f['max_d2Gdvperpdvpa_err'][:]) + max_d2Gdvpa2_err = np.copy(f['max_d2Gdvpa2_err'][:]) + max_d2Gdvperp2_err = np.copy(f['max_d2Gdvperp2_err'][:]) + expected_diff = np.copy(f['expected_diff'][:]) + expected_integral = np.copy(f['expected_integral'][:]) + + nelement_string = "N_{\\rm EL}" #\\scriptscriptstyle + ngrid_string = "N_{\\rm GR}" + + + file = filename+".plots.pdf" + pdf = PdfPages(file) + + + + marker_list = ['r--o','b--^','g--.','m--x','c--v','b','k'] + Infnorm_list = [max_dHdvpa_err, + max_dHdvperp_err,max_d2Gdvperpdvpa_err,max_d2Gdvpa2_err,max_d2Gdvperp2_err, + expected_diff, + expected_integral] + nelements = [nelement_list for item in Infnorm_list] + ylab_list = ["$\\epsilon_{\\infty}(d H / d v_{||})$","$\\epsilon_{\\infty}(d H / d v_{\\perp})$", + "$\\epsilon_{\\infty}(d^2 G / d v_{\\perp} d v_{||})$", + "$\\epsilon_{\\infty}(d^2 G / d v^2_{||})$", + "$\\epsilon_{\\infty}(d^2 G / d v^2_{\\perp})$", + "$(1/"+nelement_string+")^{"+ngrid_string+"-1}$", + "$(1/"+nelement_string+")^{"+ngrid_string+"+1}$"] + plot_1d_loglog_list_pdf (nelements,Infnorm_list,marker_list,"$"+ nelement_string+"$", pdf, + title='',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = nelement_list, yticks = None, + markersize=10, legend_title="", use_legend=True,loc_opt='lower left', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.05, 0.05), legend_fontsize=15, ncol_opt=1) + + pdf.close() + print(file) + + f.close() + return None + +workdir = "" +filename = workdir + "moment_kinetics_collisions/fkpl_integration_error_data_ngrid_5_ncore_1.h5" +get_fkpl_integration_error_data(filename) diff --git a/xarray_post_processing/plot_many_collisions.py b/xarray_post_processing/plot_many_collisions.py new file mode 100644 index 000000000..6a44cb83a --- /dev/null +++ b/xarray_post_processing/plot_many_collisions.py @@ -0,0 +1,267 @@ +# -*- coding: utf-8 -*- +import matplotlib +# Force matplotlib to not use any Xwindows backend. +matplotlib.use('Agg') # this line allows plots to be made without using a display environment variable +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np +import toml as tml +import h5py +from xarray_mk_utils import grid_data, wgts_data +from xarray_mk_utils import dynamic_data +from plot_mk_utils import plot_1d_list_pdf, plot_1d_loglog_list_pdf +from plot_mk_utils import plot_1d_semilog_list_pdf, plot_2d_pdf + + +def plot_ff_norms_with_vspace(filename,ff,ffm,vpagrid,vperpgrid): + # plot infinity norm + pdffile = filename+".ffplots.pdf" + pdf = PdfPages(pdffile) + ffplot = np.abs(ff[:,:]-ffm[:,:]) + plot_2d_pdf(vpagrid,vperpgrid,ffplot,pdf,title="$|F-F_M|$",ylab="$v_{\\perp}$",xlab="$v_{||}$") + pdf.close() + print("Saving figures: "+pdffile) + return None + +def get_time_evolving_data(filename): + print(filename) + nz_global, nz_local, zgrid = grid_data(filename,"z") + nr_global, nr_local, rgrid = grid_data(filename,"r") + nvpa_global, nvpa_local, vpagrid = grid_data(filename,"vpa") + nvperp_global, nvperp_local, vperpgrid = grid_data(filename,"vperp") + vpawgts = wgts_data(filename,"vpa") + vperpwgts = wgts_data(filename,"vperp") + time, time_present = dynamic_data(filename,"time") + ff, ff_present = dynamic_data(filename,"f") + dsdt, dsdt_present = dynamic_data(filename,"entropy_production") + density, density_present = dynamic_data(filename,"density") + parallel_flow, parallel_flow_present = dynamic_data(filename,"parallel_flow") + parallel_pressure, parallel_pressure_present = dynamic_data(filename,"parallel_pressure") + perpendicular_pressure, perpendicular_pressure_present = dynamic_data(filename,"perpendicular_pressure") + if parallel_flow_present and perpendicular_pressure_present: + pressure = (2.0*perpendicular_pressure + parallel_pressure)/3.0 + else: + pressure = None + thermal_speed, thermal_speed_present = dynamic_data(filename,"thermal_speed") + ntime = time.size + nvperp = vperpgrid.size + nvpa = vpagrid.size + #print(dsdt) + #print(density) + #print(parallel_flow) + #print(np.shape(thermal_speed)) + ffm = np.copy(ff) + for it in range(0,ntime): + for ivperp in range(0,nvperp): + for ivpa in range(0,nvpa): + vth = thermal_speed[it,0,0,0] + v2 = ((vpagrid[ivpa]-parallel_flow[it,0,0,0])/vth)**2 + (vperpgrid[ivperp]/vth)**2 + ffm[it,0,0,0,ivperp,ivpa] = density[it,0,0,0]*np.exp(-v2)/(vth**3) + + L2fm = np.copy(dsdt) + L2denom = np.sum(vperpwgts[:])*np.sum(vpawgts[:]) + #print(np.shape(L2fm)) + #print(np.shape(ff)) + for it in range(0,ntime): + L2fm[it,0,0,0] = 0.0 + for ivperp in range(0,nvperp): + for ivpa in range(0,nvpa): + L2fm[it,0,0,0] += vperpwgts[ivperp]*vpawgts[ivpa]*(ff[it,0,0,0,ivperp,ivpa]-ffm[it,0,0,0,ivperp,ivpa])**2 + #continue + L2fm[it,0,0,0] = np.sqrt(L2fm[it,0,0,0]/L2denom) + + Inffm = np.copy(dsdt) + for it in range(0,ntime): + Inffm[it,0,0,0] = 0.0 + Inffm[it,0,0,0] = np.max(np.abs(ff[it,0,0,0,:,:]-ffm[it,0,0,0,:,:])) + #plot_ff_norms_with_vspace(filename,ff[it,0,0,0,:,:],ffm[it,0,0,0,:,:],vpagrid,vperpgrid) + print("delta n: ", density[-1,0,0,0]-density[0,0,0,0]) + print("delta u: ", parallel_flow[-1,0,0,0]-parallel_flow[0,0,0,0]) + print("delta vth: ", thermal_speed[-1,0,0,0]-thermal_speed[0,0,0,0]) + #print("L2fm(t): ",L2fm[::50,0,0,0]," time: ",time[::50]) + print("L2fm: ", L2fm[0,0,0,0]," ",L2fm[-1,0,0,0]) + #print("Inffm(t): ",Inffm[::50,0,0,0]," time: ",time[::50]) + print("Inffm: ", Inffm[0,0,0,0]," ",Inffm[-1,0,0,0]) + return time, dsdt[:,0,0,0], L2fm[:,0,0,0], Inffm[:,0,0,0], density[:,0,0,0], parallel_flow[:,0,0,0], thermal_speed[:,0,0,0], pressure[:,0,0,0], vpagrid, vperpgrid, ff[-1,0,0,0,:,:], ffm[-1,0,0,0,:,:] + +def save_plot_data(filename, time, dSdt, L2norm, Infnorm, dens, upar, vth, pres, + vpagrid, vperpgrid, ff, ffm): + f = h5py.File(filename+".hdf5", "w") + f.create_dataset("time",data=time) + f.create_dataset("dSdt",data=dSdt) + f.create_dataset("L2norm",data=L2norm) + f.create_dataset("Infnorm",data=Infnorm) + f.create_dataset("dens",data=dens) + f.create_dataset("upar",data=upar) + f.create_dataset("vth",data=vth) + f.create_dataset("pres",data=pres) + f.create_dataset("vpagrid",data=vpagrid) + f.create_dataset("vperpgrid",data=vperpgrid) + f.create_dataset("ff",data=ff) + f.create_dataset("ffm",data=ffm) + f.close() + return None + +def load_plot_data(filename): + f = h5py.File(filename+".hdf5", "r") + time = np.copy(f["time"][:]) + dSdt = np.copy(f["dSdt"][:]) + L2norm = np.copy(f["L2norm"][:]) + Infnorm = np.copy(f["Infnorm"][:]) + dens = np.copy(f["dens"][:]) + upar = np.copy(f["upar"][:]) + vth = np.copy(f["vth"][:]) + pres = np.copy(f["pres"][:]) + vpagrid = np.copy(f["vpagrid"][:]) + vperpgrid = np.copy(f["vperpgrid"][:]) + ff = np.copy(f["ff"][:,:]) + ffm = np.copy(f["ffm"][:,:]) + f.close() + return time, dSdt, L2norm, Infnorm, dens, upar, vth, pres, vpagrid, vperpgrid, ff, ffm + +time_list = [] +Stime_list = [] +Mtime_list = [] +Mnoupar_time_list = [] +dSdt_list = [] +L2norm_list = [] +Infnorm_list = [] +dens_list = [] +upar_list = [] +vth_list = [] +p_list = [] +M_list = [] +Mnoupar_list = [] +M2_list = [] +M2noupar_list = [] + +#input_raw_names = ["fokker-planck-relaxation-beam-init1", +# "fokker-planck-relaxation-beam-init2", +# "fokker-planck-relaxation-beam-init3"] +#input_raw_names = ["fokker-planck-relaxation-beam-init1long", +# "fokker-planck-relaxation-beam-init2long", +# "fokker-planck-relaxation-beam-init3long"] +#input_raw_names = ["fokker-planck-relaxation-no-dfdvperp1", +# "fokker-planck-relaxation-no-dfdvperp2", +# "fokker-planck-relaxation-no-dfdvperp3"] +#input_raw_names = ["fokker-planck-relaxation-no-dfdvperp-no-conserve1", +# "fokker-planck-relaxation-no-dfdvperp-no-conserve2", +# "fokker-planck-relaxation-no-dfdvperp-no-conserve3"] +workdir = "" +input_raw_names = ["fokker-planck-relaxation-flux-bc-only1", + "fokker-planck-relaxation-flux-bc-only2", + "fokker-planck-relaxation-flux-bc-only3"] +inputname_list = [workdir+instr+".toml" for instr in input_raw_names] +outfilename_list = [workdir+instr+"/"+instr+".dfns.0.h5" for instr in input_raw_names] +process_raw_data = True + +for outfilename in outfilename_list: + savefilename = outfilename[:-10]+".processed.h5" + if process_raw_data: + time, dSdt, L2norm, Infnorm, dens, upar, vth, pres, vpagrid, vperpgrid, ff, ffm = get_time_evolving_data(outfilename) + print("Saving processed data: ",savefilename) + save_plot_data(savefilename, time, dSdt, L2norm, Infnorm, dens, upar, vth, pres, + vpagrid, vperpgrid, ff, ffm) + else: + print("Loading pre-processed data: ",savefilename) + time, dSdt, L2norm, Infnorm, dens, upar, vth, pres, vpagrid, vperpgrid, ff, ffm = load_plot_data(savefilename) + + plot_ff_norms_with_vspace(outfilename,ff,ffm,vpagrid,vperpgrid) + time_list.append(time) + Stime_list.append(time[1:]) + dSdt_list.append(dSdt[1:]) + L2norm_list.append(L2norm) + Infnorm_list.append(Infnorm) + dens_list.append(dens-dens[0]) + upar_list.append(upar-upar[0]) + vth_list.append(vth-vth[0]) + p_list.append(pres-pres[0]) + Mtime_list.append(time[:]) + Mtime_list.append(time[:]) + Mtime_list.append(time[:]) + M_list.append(np.abs(dens-dens[0])) + M_list.append(np.abs(upar-upar[0])) + M_list.append(np.abs(vth-vth[0])) + Mnoupar_time_list.append(time[:]) + Mnoupar_time_list.append(time[:]) + Mnoupar_list.append(np.abs(dens-dens[0])) + Mnoupar_list.append(np.abs(vth-vth[0])) + M2_list.append(np.abs(dens-dens[0])) + M2_list.append(np.abs(upar-upar[0])) + M2_list.append(np.abs(pres-pres[0])) + M2noupar_list.append(np.abs(dens-dens[0])) + M2noupar_list.append(np.abs(pres-pres[0])) + +file = workdir + "collisions_plots_many.pdf" +pdf = PdfPages(file) +tlabel = "$ \\nu_{ss} t $" +marker_list = ['k','r-.','b--'] +ylab_list = ["#1", "#2", "#3"] +plot_1d_semilog_list_pdf (Stime_list,dSdt_list,marker_list,tlabel, pdf, + title='$\\dot{S}$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="Resolutions", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) +plot_1d_semilog_list_pdf (time_list,L2norm_list,marker_list,tlabel, pdf, + title='$L_2(F-F_M)$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="Resolutions", use_legend=True,loc_opt='lower right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.05), legend_fontsize=15, ncol_opt=1) +plot_1d_semilog_list_pdf (time_list,Infnorm_list,marker_list,tlabel, pdf, + title='$L_{\infty}(F-F_M)$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="Resolutions", use_legend=True,loc_opt='lower right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.05), legend_fontsize=15, ncol_opt=1) + +Mylab_list = ["$|n(t)-n(0)|$ #1","$|u_{||}(t)- u_{||}(0)|$ #1","$|v_{\\rm th}(t) - v_{\\rm th}(0)|$ #1", + "$|n(t)-n(0)|$ #2","$|u_{||}(t)- u_{||}(0)|$ #2","$|v_{\\rm th}(t) - v_{\\rm th}(0)|$ #2", + "$|n(t)-n(0)|$ #3","$|u_{||}(t)- u_{||}(0)|$ #3","$|v_{\\rm th}(t) - v_{\\rm th}(0)|$ #3",] +marker_list = ['k','k-.','k--','r','r-.','r--','b','b-.','b--',] +print(len(Mtime_list)) +print(len(M_list)) +print(len(marker_list)) +print(len(Mylab_list)) +plot_1d_semilog_list_pdf (Mtime_list,M_list,marker_list,tlabel, pdf, + title='',ylab='',xlims=None,ylims=[None,10**(0)],aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = Mylab_list, + bbox_to_anchor_opt=(0.975, 0.975), legend_fontsize=15, ncol_opt=3) + +Mylab_list = ["$|n(t)-n(0)|$ #1","$|v_{\\rm th}(t) - v_{\\rm th}(0)|$ #1", + "$|n(t)-n(0)|$ #2","$|v_{\\rm th}(t) - v_{\\rm th}(0)|$ #2", + "$|n(t)-n(0)|$ #3","$|v_{\\rm th}(t) - v_{\\rm th}(0)|$ #3",] +marker_list = ['k','k--','r','r--','b','b--',] +print(len(Mnoupar_time_list)) +print(len(Mnoupar_list)) +print(len(marker_list)) +print(len(Mylab_list)) +plot_1d_semilog_list_pdf (Mnoupar_time_list,Mnoupar_list,marker_list,tlabel, pdf, + title='',ylab='',xlims=None,ylims=[None,10**(-1)],aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = Mylab_list, + bbox_to_anchor_opt=(0.975, 0.975), legend_fontsize=15, ncol_opt=3) + +Mylab_list = ["$|\\Delta n(t)|$ #1","$|\\Delta u_{||}(t)|$ #1","$|\\Delta p(t)|$ #1", + "$|\\Delta n(t)|$ #2","$|\\Delta u_{||}(t)|$ #2","$|\\Delta p(t)|$ #2", + "$|\\Delta n(t)|$ #3","$|\\Delta u_{||}(t)|$ #3","$|\\Delta p(t)|$ #3",] +marker_list = ['k','k-.','k--','r','r-.','r--','b','b-.','b--',] +print(len(Mtime_list)) +print(len(M2_list)) +print(len(marker_list)) +print(len(Mylab_list)) +#ylims = [None,10**(-7)] +#ylims = [10**(-14),10**(0)] +ylims = [None,10**(0)] +plot_1d_semilog_list_pdf (Mtime_list,M2_list,marker_list,tlabel, pdf, + title='',ylab='',xlims=None,ylims=ylims,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = Mylab_list, + bbox_to_anchor_opt=(0.825, 1.0), legend_fontsize=15, ncol_opt=3) + +Mylab_list = ["$|\\Delta n(t)|$ #1","$|\\Delta p(t)|$ #1", + "$|\\Delta n(t)|$ #2","$|\\Delta p(t)|$ #2", + "$|\\Delta n(t)|$ #3","$|\\Delta p(t)|$ #3",] +marker_list = ['k','k--','r','r--','b','b--',] +print(len(Mnoupar_time_list)) +print(len(M2noupar_list)) +print(len(marker_list)) +print(len(Mylab_list)) +plot_1d_semilog_list_pdf (Mnoupar_time_list,M2noupar_list,marker_list,tlabel, pdf, + title='',ylab='',xlims=None,ylims=[None,10**(-1)],aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = Mylab_list, + bbox_to_anchor_opt=(0.8, 0.975), legend_fontsize=15, ncol_opt=3) +pdf.close() +print("Saving figure: "+file) diff --git a/xarray_post_processing/plot_mk_utils.py b/xarray_post_processing/plot_mk_utils.py new file mode 100644 index 000000000..01d52ecb9 --- /dev/null +++ b/xarray_post_processing/plot_mk_utils.py @@ -0,0 +1,148 @@ +import matplotlib.pyplot as plt +from matplotlib import rcParams +import numpy as np +plt.rc('text', usetex=False) +plt.rc('font', family='serif') +plt.rc('font', size=20) +rcParams.update({'text.latex.preamble' : r'\usepackage{bm}'}) +rcParams.update({'figure.autolayout': True}) + +def plot_1d_list_pdf (xlist,ylist,marker_list,xlab, pdf, + title='',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, xticks_labels=None, yticks = None, + markersize=5, legend_title="", use_legend=False,loc_opt='upper right', ylab_list = None, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=10, ncol_opt=1, + legend_shadow=False,legend_frame=False, vlines = None,hlines = None,marker_fill_style = None, + cartoon=False, linewidth=None, texts = None, slines=None): + + fig=plt.figure(figsize=(aspx,aspy)) + nlist = len(ylist) + if(ylab_list is None): + ylab_list = [None for i in range(0,nlist)] + for iy in range(0,nlist): + plt.plot(xlist[iy],ylist[iy],marker_list[iy],markersize=markersize,label=ylab_list[iy], + fillstyle = marker_fill_style, linewidth = linewidth) + plt.xlabel(xlab) + if len(ylab) > 0: + plt.ylabel(ylab) + if len(title) > 0: + plt.title(title) + if(not xlims is None): + plt.xlim(xlims[0],xlims[1]) + if(not ylims is None): + plt.ylim(ylims[0],ylims[1]) + if(not vlines is None): + for xin,xlabel,xcolor,xlinestyle in vlines: + plt.axvline(x=xin, label=xlabel, color=xcolor,linestyle=xlinestyle,linewidth=linewidth) + if(not hlines is None): + for yin,ylabel,ycolor,ylinestyle in hlines: + plt.axhline(y=yin, label=ylabel, color=ycolor,linestyle=ylinestyle,linewidth=linewidth) + if(not texts is None): + for xin, yin, textin in texts: + print(xin,yin,textin) + plt.text(xin,yin,textin) + if (not slines is None): + for m,c,marker,label in slines: + plt.plot(xlist[0],m*xlist[0]+c,marker,label=label) + if(use_legend): + plt.legend(title=legend_title,loc=loc_opt, bbox_to_anchor=bbox_to_anchor_opt, + fontsize=legend_fontsize, frameon=legend_frame, handlelength=1, labelspacing=0.5, + ncol=ncol_opt, columnspacing = 0.5 , handletextpad = 0.5, shadow=legend_shadow) + if(not xticks is None): + plt.xticks(xticks) + if(not yticks is None): + plt.yticks(yticks) + if (cartoon): + plt.tick_params(top='off', bottom='off', left='off', right='off', labelleft='off', labelbottom='off') + plt.box(False) + pdf.savefig(fig)# pdf is the object of the current open PDF file to which the figures are appended + plt.close (fig) + return + +def plot_1d_semilog_list_pdf (xlist,ylist,marker_list,xlab, pdf, + title='',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=False,loc_opt='upper right', ylab_list = None, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=10, ncol_opt=1, + legend_shadow=False,legend_frame=False): + + fig=plt.figure(figsize=(aspx,aspy)) + nlist = len(ylist) + if(ylab_list is None): + ylab_list = [None for i in range(0,nlist)] + for iy in range(0,nlist): + plt.semilogy(xlist[iy],ylist[iy],marker_list[iy],markersize=markersize,label=ylab_list[iy]) + plt.xlabel(xlab) + if len(ylab) > 0: + plt.ylabel(ylab) + if len(title) > 0: + plt.title(title) + if(not xlims is None): + plt.xlim(xlims[0],xlims[1]) + if(not ylims is None): + plt.ylim(ylims[0],ylims[1]) + if(use_legend): + plt.legend(title=legend_title,loc=loc_opt, bbox_to_anchor=bbox_to_anchor_opt, + fontsize=legend_fontsize, frameon=legend_frame, handlelength=1, labelspacing=0.5, + ncol=ncol_opt, columnspacing = 0.5 , handletextpad = 0.5, shadow=legend_shadow) + if(not xticks is None): + plt.xticks(xticks) + if(not yticks is None): + plt.yticks(yticks) + pdf.savefig(fig)# pdf is the object of the current open PDF file to which the figures are appended + plt.close (fig) + return + +def plot_1d_loglog_list_pdf (xlist,ylist,marker_list,xlab, pdf, + title='',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=False,loc_opt='upper right', ylab_list = None, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=10, ncol_opt=1, + legend_shadow=False,legend_frame=False): + + fig=plt.figure(figsize=(aspx,aspy)) + nlist = len(ylist) + if(ylab_list is None): + ylab_list = [None for i in range(0,nlist)] + for iy in range(0,nlist): + plt.loglog(xlist[iy],ylist[iy],marker_list[iy],markersize=markersize,label=ylab_list[iy]) + plt.xlabel(xlab) + if len(ylab) > 0: + plt.ylabel(ylab) + if len(title) > 0: + plt.title(title) + if(not xlims is None): + plt.xlim(xlims[0],xlims[1]) + if(not ylims is None): + plt.ylim(ylims[0],ylims[1]) + if(use_legend): + plt.legend(title=legend_title,loc=loc_opt, bbox_to_anchor=bbox_to_anchor_opt, + fontsize=legend_fontsize, frameon=legend_frame, handlelength=1, labelspacing=0.5, + ncol=ncol_opt, columnspacing = 0.5 , handletextpad = 0.5, shadow=legend_shadow) + if(not xticks is None): + plt.xticks([]) + plt.minorticks_off() + #print(plt.xticks()) + plt.xticks(xticks,[str(tick) for tick in xticks]) + #print(plt.xticks()) + + if(not yticks is None): + plt.yticks(yticks) + pdf.savefig(fig)# pdf is the object of the current open PDF file to which the figures are appended + plt.close (fig) + return + +def plot_2d_pdf(x,y,z,pdf,title="",ylab="",xlab=""): + + # make data + X, Y = np.meshgrid(x, y) + levels = np.linspace(z.min(), z.max(), 7) + + # plot + fig = plt.figure() + + plt.contourf(X, Y, z, levels=levels) + plt.colorbar() + plt.title(title) + plt.xlabel(xlab) + plt.ylabel(ylab) + pdf.savefig(fig) + plt.close(fig) + return None diff --git a/xarray_post_processing/plot_sd.py b/xarray_post_processing/plot_sd.py new file mode 100644 index 000000000..dfc9548a1 --- /dev/null +++ b/xarray_post_processing/plot_sd.py @@ -0,0 +1,176 @@ +# -*- coding: utf-8 -*- +import matplotlib +# Force matplotlib to not use any Xwindows backend. +matplotlib.use('Agg') # this line allows plots to be made without using a display environment variable +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np +import toml as tml +from xarray_mk_utils import grid_data +from xarray_mk_utils import dynamic_data +from plot_mk_utils import plot_1d_list_pdf, plot_1d_loglog_list_pdf + +def get_sd_plot_data(filename): + # return grid, moments, and pdf at the last timestep of the simulation + # where the pdf is written + # assume that there is single output file, from a simulation + # using parallel HDF5 or a single shared-memory region + + nz_global, nz_local, zgrid = grid_data(filename,"z") + nr_global, nr_local, rgrid = grid_data(filename,"r") + nvpa_global, nvpa_local, vpagrid = grid_data(filename,"vpa") + nvperp_global, nvperp_local, vperpgrid = grid_data(filename,"vperp") + + time, time_present = dynamic_data(filename,"time") + ff, ff_present = dynamic_data(filename,"f") + + ntime = time.size + + #print("z: ",zgrid) + #print("r: ",rgrid) + #print("vpa: ",vpagrid) + #print("vperp: ",vperpgrid) + + # return the data of interest only + it = ntime - 1 + ivpa = nvpa_global//2+1 + ivperp = 0 + ir = 0 + iz = 0 # lower wall plate + ispec = 0 # single species + + return vperpgrid, ff[it,ispec,ir,iz,:,ivpa], vpagrid, ff[it,ispec,ir,iz,ivperp,:] + +def get_sd_input_data(filename): + print(filename) + with open(filename, 'r') as file: + inputdata = tml.load(file) + key = "fokker_planck_collisions" + print(inputdata[key]) + ni = inputdata[key]["sd_density"] + Ti = inputdata[key]["sd_temp"] + Te = Ti + # only a single charge ion is evolved, hence + # we interpret species i as alpha here + Zalpha = inputdata[key]["Zi"] + # Zi must be absolute / relative to proton charge + # whereas the fixed Maxwellian ion charge number is given here + Zi = inputdata[key]["sd_q"] + ne = Zi*ni + 1.0*Zalpha # quasineutrality to determine electron density + # initial alpha density must be unity + mi = inputdata[key]["sd_mi"] + me = inputdata[key]["sd_me"] + nuref = inputdata[key]["nuii"] + # compute critical speed vc3/cref^3, cref = sqrt(2 Tref/mref) -> factor of 1/ 2 sqrt 2 + vc3 = (np.sqrt(2.0)/4.0)*3.0*np.sqrt(np.pi/2.0)*(Zi**2)*((Te)**1.5)*(ni/ne)/(np.sqrt(me)*mi) + + key = "ion_source" + print(inputdata[key]) + v0 = inputdata[key]["source_v0"] + Salpha = inputdata[key]["source_strength"] + # use that nuref = gamma_alphaalpha nref / mref^2 cref^3, with cref = sqrt(2Tref/mref) and alphas the reference species + # gamma_alphaalpha = 2 pi Zalpha^4 e^4 ln Lambda/ (4pi epsilon0)^2 + # and nu_alphae = (4/(3sqrt(2pi))) gamma_alphae ne Te^(-3/2) sqrt(me)/m_alpha + nualphae = nuref*(8.0/3.0)*(1.0/np.sqrt(np.pi))*ne*np.sqrt(me)*(Te**(-1.5))*(Zalpha**(2)) + amplitude = (np.sqrt(np.pi)/4.0)*Salpha/nualphae # pi^3/2 * (1/4 pi) factor had pi^3/2 due to normalisation of integration and pdf + return v0, vc3, amplitude + +workdir = "" +input_filename_list = [workdir+"/excalibur/moment_kinetics_gyro/runs/fokker-planck-relaxation-example-4.toml"] +filename_dfns_list = [workdir+"/excalibur/moment_kinetics_gyro/runs/fokker-planck-relaxation-example-4/fokker-planck-relaxation-example-4.dfns.0.h5",] +identity = "example-4" + +vpagrid_list = [] +vperpgrid_list = [] +ff_list = [] +logff_list = [] +ffvpa_list = [] +logffvpa_list = [] + +for ifile, filename_dfn in enumerate(filename_dfns_list): + + v0, vc3, amplitude = get_sd_input_data(input_filename_list[ifile]) + + vperpgrid, ff, vpagrid, ffvpa = get_sd_plot_data(filename_dfn) + vperpgrid_list.append(vperpgrid) + ff_list.append(ff) + logff_list.append(np.log(np.abs(ff)+1.0e-15)) + vpagrid_list.append(vpagrid) + ffvpa_list.append(ffvpa) + logffvpa_list.append(np.log(np.abs(ffvpa)+1.0e-15)) + + # compute a slowing down distribution for comparison from an analytical formula + + ff_sd = np.copy(vperpgrid) + nvperp = vperpgrid.size + vc3test = 3.0*np.sqrt(np.pi/2.0)*((0.01)**1.5)*(1.0/(np.sqrt(2.7e-4)*0.5))*(0.5*0.5*1.0/2.0) + #print(vc3test," ",vc3test**(1.0/3.0)) + print("vc3: ", vc3," vc: ",vc3**(1.0/3.0)) + for ivperp in range(0,nvperp): + vperp = vperpgrid[ivperp] + if vperp < v0: + ff_sd[ivperp] = 1.0/(vc3 + vperp**3.0) + #print(ivperp) + else: + ff_sd[ivperp] = 0.0 + # pick a point to normalise by + ivperp = 32 #nvperp//3 + 1 + amplitude_test=ff[ivperp]/ff_sd[ivperp] + print(amplitude_test," ",amplitude, " ", amplitude/amplitude_test) + ff_sd = ff_sd*amplitude +# ff_sd = ff_sd*amplitude_test + + vperpgrid_list.append(vperpgrid) + ff_list.append(ff_sd) + logff_list.append(np.log(np.abs(ff_sd)+1.0e-15)) + + ff_sd_vpa = np.copy(vpagrid) + nvpa = vpagrid.size + for ivpa in range(0,nvpa): + vpa = vpagrid[ivpa] + if np.abs(vpa) < v0: + ff_sd_vpa[ivpa] = 1.0/(vc3 + np.abs(vpa)**3.0) + #print(ivpa) + else: + ff_sd_vpa[ivpa] = 0.0 + # pick a point to normalise by + ivpa = 96# nvpa//2 + nvpa//6 + 1 + amplitude_test=ffvpa[ivperp]/ff_sd_vpa[ivperp] + print(amplitude_test," ",amplitude, " ", amplitude/amplitude_test) + ff_sd_vpa = ff_sd_vpa*amplitude +# ff_sd_vpa = ff_sd_vpa*amplitude_test + + vpagrid_list.append(vpagrid) + ffvpa_list.append(ff_sd_vpa) + logffvpa_list.append(np.log(np.abs(ff_sd_vpa)+1.0e-15)) + +marker_list = ['k','b','r','g','c','y'] +ylab_list = ["Num","SD"] +file = workdir + "excalibur/moment_kinetics_gyro/sd_scan_"+str(identity)+".pdf" +pdf = PdfPages(file) + +# plot ff +plot_1d_list_pdf (vperpgrid_list,ff_list,marker_list,"$v_{\\perp}$", pdf, + title='$f(v_{\\|}=0,v_{\\perp})$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) +# plot logff +plot_1d_list_pdf (vperpgrid_list,logff_list,marker_list,"$v_{\\perp}$", pdf, + title='$\\ln|f(v_{\\|}=0,v_{\\perp})|$',ylab='',xlims=None,ylims=[-10.0,1.1*np.max(logff_list[0])],aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) + +# plot ffvpa +#print(vpagrid_list,ffvpa_list) +plot_1d_list_pdf (vpagrid_list,ffvpa_list,marker_list,"$v_{\\|}$", pdf, + title='$f(v_{\\|},v_{\\perp}=0)$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) +# plot logffvpa +plot_1d_list_pdf (vpagrid_list,logffvpa_list,marker_list,"$v_{\\|}$", pdf, + title='$\\ln|f(v_{\\|},v_{\\perp}=0)|$',ylab='',xlims=None,ylims=[-10.0,1.1*np.max(logffvpa_list[0])],aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) + + +pdf.close() +print("Saving figure: "+file) diff --git a/xarray_post_processing/plot_wall.py b/xarray_post_processing/plot_wall.py new file mode 100644 index 000000000..5eb3c5865 --- /dev/null +++ b/xarray_post_processing/plot_wall.py @@ -0,0 +1,150 @@ +# -*- coding: utf-8 -*- +import matplotlib +# Force matplotlib to not use any Xwindows backend. +matplotlib.use('Agg') # this line allows plots to be made without using a display environment variable +from matplotlib.backends.backend_pdf import PdfPages +import numpy as np + +from xarray_mk_utils import grid_data +from xarray_mk_utils import dynamic_data +from plot_mk_utils import plot_1d_list_pdf, plot_1d_loglog_list_pdf + +def get_wall_plot_data(filename): + # return grid, moments, and pdf at the last timestep of the simulation + # where the pdf is written + # assume that there is single output file, from a simulation + # using parallel HDF5 or a single shared-memory region + + nz_global, nz_local, zgrid = grid_data(filename,"z") + nr_global, nr_local, rgrid = grid_data(filename,"r") + nvpa_global, nvpa_local, vpagrid = grid_data(filename,"vpa") + nvperp_global, nvperp_local, vperpgrid = grid_data(filename,"vperp") + + time, time_present = dynamic_data(filename,"time") + ff, ff_present = dynamic_data(filename,"f") + Ez, Ez_present = dynamic_data(filename,"Ez") + phi, phi_present = dynamic_data(filename,"phi") + cil, cil_present = dynamic_data(filename,"chodura_integral_lower") + + ntime = time.size + + #print("z: ",zgrid) + #print("r: ",rgrid) + #print("vpa: ",vpagrid) + #print("vperp: ",vperpgrid) + + # return the data of interest only + it = ntime - 1 + ivperp = 0 + ir = 0 + iz = 0 # lower wall plate + ispec = 0 # single species + return zgrid, vpagrid, ff[it,ispec,ir,iz,ivperp,:], Ez[it,ir,:], phi[it,ir,:], cil[it,ir] + +workdir = "" +filename = workdir + "moment_kinetics_newgeo/runs/wall-bc_cheb/wall-bc_cheb.moments.0.h5" +filename_dfns_list = ["moment_kinetics_newgeo/runs/wall-bc_cheb_epsz1/wall-bc_cheb_epsz1.dfns.0.h5", + "moment_kinetics_newgeo/runs/wall-bc_cheb_epsz0.1/wall-bc_cheb_epsz0.1.dfns.0.h5", + "moment_kinetics_newgeo/runs/wall-bc_cheb_epsz0.01/wall-bc_cheb_epsz0.01.dfns.0.h5", + "moment_kinetics_newgeo/runs/wall-bc_cheb_epsz0.001/wall-bc_cheb_epsz0.001.dfns.0.h5", + "moment_kinetics_newgeo/runs/wall-bc_cheb_epsz0/wall-bc_cheb_epsz0.dfns.0.h5"] + +zgrid_list = [] +vpagrid_list = [] +ff_list = [] +ff_over_vpa2_list = [] +Ez_list = [] +phi_list = [] +cil_list = [] +logphi_list = [] +logEz_list = [] +logz_list = [] +for filename_dfn in filename_dfns_list: + zgrid, vpagrid, ff, Ez, phi, cil = get_wall_plot_data(workdir+filename_dfn) + nz = zgrid.size + zgrid_list.append(zgrid) + vpagrid_list.append(vpagrid) + ff_list.append(ff) + Ez_list.append(Ez) + phi_list.append(phi) + cil_list.append(cil) + nzlog = nz//12 + #print(zgrid[-1-nzlog:-1]) + logz_list.append(np.log(0.5-zgrid[-1-nzlog:-1])) + logphi_list.append(np.log(phi[-1-nzlog:-1]-phi[-1])) + logEz_list.append(np.log(Ez[-1-nzlog:-1])) + nvpa = vpagrid.size + vpafunc = np.zeros(nvpa) + deltavpa = np.amin(vpagrid[1:nvpa]-vpagrid[:nvpa-1]) + #print(deltavpa) + for ivpa in range(0,nvpa): + if np.abs(vpagrid[ivpa]) > 0.5*deltavpa: + vpafunc[ivpa] = 1.0/(vpagrid[ivpa]**2) + ff_over_vpa2_list.append(ff*vpafunc) + +#print(logz_list) +#print(logEz_list) +#print(logphi_list) +#print(ff_over_vpa2_list) +epsz_values = [1.0,0.1,0.01,0.001,0.0] +marker_list = ['k','b','r','g','c','y'] +ylab_list = [str(epsz) for epsz in epsz_values] +file = workdir + "moment_kinetics_newgeo/wall_boundary_cutoff_scan.pdf" +pdf = PdfPages(file) + +# plot ff +plot_1d_list_pdf (vpagrid_list,ff_list,marker_list,"$v_{\\|\\|}$", pdf, + title='$f(z_{\\rm wall-},v_{\\|\\|})$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="$\\epsilon_z$", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) +# plot ff/vpa2 +plot_1d_list_pdf (vpagrid_list,ff_over_vpa2_list,marker_list,"$v_{\\|\\|}$", pdf, + title='$f(z_{\\rm wall-},v_{\\|\\|})/v_{\\|\\|}^2$',ylab='',xlims=None,ylims=[-0.1,5.0],aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="$\\epsilon_z$", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1) +# plot Bohm condition +ylist = [np.array(cil_list)] +xlist = [np.array(epsz_values)] +plot_1d_list_pdf (xlist,ylist,["kx--"],"$\\epsilon_z$", pdf, + title='$(T_{\\rm e}/2 n) \\int (f/v_{\\|\\|}^2) d v_{\\|\\|}/\\sqrt{\\pi} $',ylab='', + xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=15, legend_title="", use_legend=False,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.25, 1.0), legend_fontsize=15, ncol_opt=1, hlines = [[1.0,"","r","--"]]) +# plot phi +plot_1d_list_pdf (zgrid_list,phi_list,marker_list,"$z/L_z$", pdf, + title='$\\phi(z)$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="$\\epsilon_z$", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.55, 0.85), legend_fontsize=15, ncol_opt=1) +# plot Ez +plot_1d_list_pdf (zgrid_list,Ez_list,marker_list,"$z/L_z$", pdf, + title='$E_z(z)$',ylab='',xlims=None,ylims=None,aspx=9,aspy=6, xticks = None, yticks = None, + markersize=5, legend_title="$\\epsilon_z$", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.25, 1.0), legend_fontsize=15, ncol_opt=1) +# plot log phi +plot_1d_list_pdf (logz_list,logphi_list,marker_list,"$\\ln(0.5 - z/L_z)$", pdf, + title='$\\ln (\\phi(z)-\\phi_{\\rm wall})$',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1, + legend_shadow=False,legend_frame=False,slines=[[0.5,0.0,'k--'," 0.5 log dz"],[0.6666, 1.5, 'b--',"2 log dz/3 + 1.5"]]) +# plot log Ez +plot_1d_list_pdf (logz_list,logEz_list,marker_list,"$\\ln(0.5 - z/L_z)$", pdf, + title='$\\ln E_z(z)$',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, yticks = None, + markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, + bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1, + legend_shadow=False,legend_frame=False,slines=[[-0.5,0.0,'k--'," -0.5 log dz"],[-0.333, 0.9,'b--',"0.9 - log dz/3"]]) + +# plot log phi +#plot_1d_loglog_list_pdf (logz_list,logphi_list,marker_list,"$0.5 - z/L_z$", pdf, +# title='$\\phi(z)-\\phi_{\\rm wall}$',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, yticks = None, +# markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, +# bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1, +# legend_shadow=False,legend_frame=False) +# plot log Ez +#plot_1d_loglog_list_pdf (logz_list,logEz_list,marker_list,"$0.5 - z/L_z$", pdf, +# title='$E_z(z)$',ylab='',xlims=None,ylims=None,aspx=12,aspy=8, xticks = None, yticks = None, +# markersize=5, legend_title="", use_legend=True,loc_opt='upper right', ylab_list = ylab_list, +# bbox_to_anchor_opt=(0.95, 0.95), legend_fontsize=15, ncol_opt=1, +# legend_shadow=False,legend_frame=False) + +pdf.close() +print("Saving figure: "+file) diff --git a/xarray_post_processing/xarray_mk_utils.py b/xarray_post_processing/xarray_mk_utils.py new file mode 100644 index 000000000..c5106e197 --- /dev/null +++ b/xarray_post_processing/xarray_mk_utils.py @@ -0,0 +1,33 @@ +import xarray as xr + + +def read_variable(dataset,varstring): + try: + var=dataset[varstring].data + var_present=True + except KeyError: + print('INFO: '+varstring+' not found in data file') + var=None + var_present = False + return var, var_present + +def grid_data(filename,coord): + dataset = xr.open_dataset(filename,group ="coords/"+coord) + n_global = dataset["n_global"].data + n_local = dataset["n_local"].data + grid = dataset["grid"].data + dataset.close() + return n_global, n_local, grid + +def wgts_data(filename,coord): + dataset = xr.open_dataset(filename,group ="coords/"+coord) + wgts = dataset["wgts"].data + dataset.close() + return wgts + +def dynamic_data(filename,varstring): + dataset = xr.open_dataset(filename,group ="dynamic_data/") + var, var_present = read_variable(dataset,varstring) + return var, var_present + + From d40a160e5c83cbd088c3ece2bd15dd637f12ae84 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Tue, 6 Aug 2024 13:03:39 +0100 Subject: [PATCH 2/3] Move xarray_post_processing files and add README.md and requirements.txt. --- .../xarray_post_processing/README.md | 19 +++++++++++++++++++ .../plot_error_data.py | 0 .../plot_integration_error_data.py | 0 .../plot_many_collisions.py | 0 .../xarray_post_processing}/plot_mk_utils.py | 0 .../xarray_post_processing}/plot_sd.py | 0 .../xarray_post_processing}/plot_wall.py | 0 .../xarray_post_processing/requirements.txt | 5 +++++ .../xarray_mk_utils.py | 0 9 files changed, 24 insertions(+) create mode 100644 publication_inputs/xarray_post_processing/README.md rename {xarray_post_processing => publication_inputs/xarray_post_processing}/plot_error_data.py (100%) rename {xarray_post_processing => publication_inputs/xarray_post_processing}/plot_integration_error_data.py (100%) rename {xarray_post_processing => publication_inputs/xarray_post_processing}/plot_many_collisions.py (100%) rename {xarray_post_processing => publication_inputs/xarray_post_processing}/plot_mk_utils.py (100%) rename {xarray_post_processing => publication_inputs/xarray_post_processing}/plot_sd.py (100%) rename {xarray_post_processing => publication_inputs/xarray_post_processing}/plot_wall.py (100%) create mode 100644 publication_inputs/xarray_post_processing/requirements.txt rename {xarray_post_processing => publication_inputs/xarray_post_processing}/xarray_mk_utils.py (100%) diff --git a/publication_inputs/xarray_post_processing/README.md b/publication_inputs/xarray_post_processing/README.md new file mode 100644 index 000000000..3275baee0 --- /dev/null +++ b/publication_inputs/xarray_post_processing/README.md @@ -0,0 +1,19 @@ + +## xarray & h5py plotting scripts for publication quality figures + +This directory contains python scripts for making publication quality figures. +We briefly describe the contents of the files. + +* `plot_mk_utils.py`: A series of plotting functions using matplotlib and pyplot. + +* `xarray_mk_utils.py`: A series of utility functions for reading data from `moment_kinetics` output files. + +* `plot_wall.py`: A script for comparing sheath-boundary simulations. + +* `plot_sd.py`: A script for comparing the numerical solution from `moment_kinetics` to the analytical slowing-down solution. + +* `plot_many_collisions.py`: A script for comparing simulations of the relaxation to the Maxwellian distribution in the presence of self collisions. + +* `plot_error_data.py` and `plot_integration_error_data.py`: Scripts for plotting data produced by the evaluation tests of the Fokker-Planck collision operator. + +The `requirements.txt` file provides a list of required modules at the last used version. diff --git a/xarray_post_processing/plot_error_data.py b/publication_inputs/xarray_post_processing/plot_error_data.py similarity index 100% rename from xarray_post_processing/plot_error_data.py rename to publication_inputs/xarray_post_processing/plot_error_data.py diff --git a/xarray_post_processing/plot_integration_error_data.py b/publication_inputs/xarray_post_processing/plot_integration_error_data.py similarity index 100% rename from xarray_post_processing/plot_integration_error_data.py rename to publication_inputs/xarray_post_processing/plot_integration_error_data.py diff --git a/xarray_post_processing/plot_many_collisions.py b/publication_inputs/xarray_post_processing/plot_many_collisions.py similarity index 100% rename from xarray_post_processing/plot_many_collisions.py rename to publication_inputs/xarray_post_processing/plot_many_collisions.py diff --git a/xarray_post_processing/plot_mk_utils.py b/publication_inputs/xarray_post_processing/plot_mk_utils.py similarity index 100% rename from xarray_post_processing/plot_mk_utils.py rename to publication_inputs/xarray_post_processing/plot_mk_utils.py diff --git a/xarray_post_processing/plot_sd.py b/publication_inputs/xarray_post_processing/plot_sd.py similarity index 100% rename from xarray_post_processing/plot_sd.py rename to publication_inputs/xarray_post_processing/plot_sd.py diff --git a/xarray_post_processing/plot_wall.py b/publication_inputs/xarray_post_processing/plot_wall.py similarity index 100% rename from xarray_post_processing/plot_wall.py rename to publication_inputs/xarray_post_processing/plot_wall.py diff --git a/publication_inputs/xarray_post_processing/requirements.txt b/publication_inputs/xarray_post_processing/requirements.txt new file mode 100644 index 000000000..bae625561 --- /dev/null +++ b/publication_inputs/xarray_post_processing/requirements.txt @@ -0,0 +1,5 @@ +h5py==3.11.0 +matplotlib==3.7.3 +numpy==1.24.4 +toml==0.10.2 +xarray==2023.1.0 diff --git a/xarray_post_processing/xarray_mk_utils.py b/publication_inputs/xarray_post_processing/xarray_mk_utils.py similarity index 100% rename from xarray_post_processing/xarray_mk_utils.py rename to publication_inputs/xarray_post_processing/xarray_mk_utils.py From 1d680ce4b75924da143c89fa22830e24469bcfc8 Mon Sep 17 00:00:00 2001 From: Michael Hardman <29800382+mrhardman@users.noreply.github.com> Date: Mon, 2 Sep 2024 12:13:28 +0100 Subject: [PATCH 3/3] Change plot labels epsilon_L_2 -> epsilon_2. --- .../xarray_post_processing/plot_error_data.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/publication_inputs/xarray_post_processing/plot_error_data.py b/publication_inputs/xarray_post_processing/plot_error_data.py index 1dfac40f9..cc185de28 100644 --- a/publication_inputs/xarray_post_processing/plot_error_data.py +++ b/publication_inputs/xarray_post_processing/plot_error_data.py @@ -54,7 +54,7 @@ def get_fkpl_error_data(filename): marker_list = ['r--o','b--s','g--.','m--x','c--v','k'] C_list = [max_C_err,L2_C_err,n_err,u_err,p_err,expected_diff] nelements = [nelement_list for item in C_list] - ylab_list = ["$\\epsilon_{\\infty}(C[F,F])$", "$\\epsilon_{L_2}(C[F,F])$", + ylab_list = ["$\\epsilon_{\\infty}(C[F,F])$", "$\\epsilon_{2}(C[F,F])$", "$|\\Delta n|$", "$|\\Delta u_{||}|$", "$|\\Delta p |$", @@ -98,10 +98,10 @@ def get_fkpl_error_data(filename): #expected_diff, expected_integral] nelements = [nelement_list for item in L2norm_list] - ylab_list = ["$\\epsilon_{L_2}(d H / d v_{||})$","$\\epsilon_{L_2}(d H / d v_{\\perp})$", - "$\\epsilon_{L_2}(d^2 G / d v_{\\perp} d v_{||})$", - "$\\epsilon_{L_2}(d^2 G / d v^2_{||})$", - "$\\epsilon_{L_2}(d^2 G / d v^2_{\\perp})$", + ylab_list = ["$\\epsilon_{2}(d H / d v_{||})$","$\\epsilon_{2}(d H / d v_{\\perp})$", + "$\\epsilon_{2}(d^2 G / d v_{\\perp} d v_{||})$", + "$\\epsilon_{2}(d^2 G / d v^2_{||})$", + "$\\epsilon_{2}(d^2 G / d v^2_{\\perp})$", #"$(1/"+nelement_string+")^{"+ngrid_string+"-1}$", "$(1/"+nelement_string+")^{"+ngrid_string+"+1}$"] plot_1d_loglog_list_pdf (nelements,L2norm_list,marker_list,"$"+ nelement_string+"$", pdf,