diff --git a/README.md b/README.md index 601ca86..e8b74bc 100644 --- a/README.md +++ b/README.md @@ -12,7 +12,7 @@ For more information on the [libcloudph++](http://libcloudphxx.igf.fuw.edu.pl/) [project repository](https://github.com/igfuw/libcloudphxx) and the [project documentation](http://www.geosci-model-dev.net/8/1677/2015/). -The parcel model is written in Python 2.7. +The parcel model was written in Python 2.7 and now works with Python 3.# . # installation diff --git a/long_test/test_timestep.py b/long_test/test_timestep.py index 1087972..da84926 100644 --- a/long_test/test_timestep.py +++ b/long_test/test_timestep.py @@ -1,5 +1,6 @@ import sys sys.path.insert(0, "../") +#sys.path.insert(0, "path to /parcel/plots/comparison") sys.path.insert(0, "./") sys.path.insert(0, "plots/comparison/") @@ -15,33 +16,35 @@ import pdb """ -This set of tests checks how the timestep of the simulation affects -the maximum supersaturation (RH) and the concentration of the activated -particles (N). -The expected result is to see nearly constant RH and N for small timesteps -(.001 - .03 for this setup) and then decrease of RH and increase of N -for bigger timesteps. +This set of tests checks how the timestep of the simulation affects +the maximum supersaturation (RH) and the concentration of the activated +particles (N). +The expected result is to see nearly constant RH and N for small timesteps +(.001 - .03 for this setup) and then decrease of RH and increase of N +for bigger timesteps. """ Dt_list = [1e-3, 2e-3, 4e-3, 8e-3, 1e-2, 2e-2, 4e-2, 1e-1, 2e-1, 1.] -# runs all simulations +# runs all simulations # returns data with values of RH_max and N at the end of simulations # keep the netcdf alive for all tests @pytest.fixture(scope="module") def data(request): - # initial values + # initial values RH_init = .99999 T_init = 280. p_init = 100000. r_init = common.eps * RH_init * common.p_vs(T_init) / (p_init - RH_init * common.p_vs(T_init)) - # lists to store RH_max and N at the end of the simulation from each test run + # lists to store RH_max and N at the end of the simulation from each test run RH_list = [] N_list = [] for dt in Dt_list: + print("\nt time step", dt) outfile_nc = "timesteptest_dt=" + str(dt) + ".nc" + parcel(dt=dt, outfreq = int(100/dt), outfile = outfile_nc,\ w = 1., T_0 = T_init, p_0 = p_init, r_0 = r_init, z_max = 200, \ sd_conc = 1000, \ @@ -52,9 +55,9 @@ def data(request): f_out = netcdf.netcdf_file(outfile_nc, "r") RH_max = f_out.RH_max N_end = f_out.variables["radii_m0"][-1,0] # concentration of drops > 1e-6 m - - RH_list.append((RH_max - 1)*100) # [%] - N_list.append(N_end / 1e6) # [1/mg] + + RH_list.append((RH_max - 1)*100) # [%] + N_list.append(N_end / 1e6) # [1/mg] data = {"RH" : RH_list, "N" : N_list, "dt" : Dt_list} @@ -68,10 +71,10 @@ def removing_files(): def test_timestep_eps(data, eps=0.01, dt_lim=0.01): """ - checking if the results obtained from simulations with different timesteps - do not differ from the referential one (the one with the smallest timestep) - more than eps times - (Unitill we think of a better convergence test, the check is done for the + checking if the results obtained from simulations with different timesteps + do not differ from the referential one (the one with the smallest timestep) + more than eps times + (Unitill we think of a better convergence test, the check is done for the smallest timesteps. This is done in order to avoid too big epsilon. """ for var, val in data.items(): @@ -79,24 +82,24 @@ def test_timestep_eps(data, eps=0.01, dt_lim=0.01): if var in ["RH", "N"]: # for simulations with small timesteps for idx in range(len(data["dt"])): - if data["dt"][idx] < dt_lim: + if data["dt"][idx] < dt_lim: # assert that the results are close to the one with the smallest timestep - assert np.isclose(val[idx], val[0], atol=0, rtol=eps), str(val[idx]) + str(val[0]) + assert np.isclose(val[idx], val[0], atol=0, rtol=eps), str(val[idx]) + str(val[0]) @pytest.mark.parametrize("dt", Dt_list) def test_timestep_diff(data, dt, eps=2.e-4): """ - checking if the results are close to the referential ones + checking if the results are close to the referential ones (stored in refdata folder) """ filename = "timesteptest_dt=" + str(dt) + ".nc" f_test = netcdf.netcdf_file(filename, "r") - f_ref = netcdf.netcdf_file(os.path.join("long_test/refdata", filename), "r") + f_ref = netcdf.netcdf_file(os.path.join("/home/piotr/Piotr/IGF/parcel/long_test/refdata", filename), "r") for var in ["t", "z", "th_d", "T", "p", "r_v", "rhod"]: assert np.isclose(f_test.variables[var][:], f_ref.variables[var][:], atol=0, rtol=eps).all(), "differs e.g. " + str(var) + "; max(ref diff) = " + str(np.where(f_ref.variables[var][:] != 0., abs((f_test.variables[var][:]-f_ref.variables[var][:])/f_ref.variables[var][:]), 0.).max()) - + def test_timestep_plot(data): timestep_plot(data, output_folder="plots/outputs") diff --git a/parcel.py b/parcel.py index f274658..523c388 100755 --- a/parcel.py +++ b/parcel.py @@ -5,6 +5,7 @@ #sys.path.insert(0, "../libcloudphxx/build/bindings/python/") #sys.path.insert(0, "../../../libcloudphxx/build/bindings/python/") #sys.path.insert(0, "/usr/local/lib/site-python/") +#sys.path.insert(0, "path_to/lib/python3/dist-packages") # TEMP TODO TEMP TODO !!! from argparse import ArgumentParser, RawTextHelpFormatter @@ -155,17 +156,19 @@ def _output_bins(fout, t, micro, opts, spectra): for dim, dct in spectra.items(): for bin in range(dct["nbin"]): if dct["drwt"] == 'wet': - micro.diag_wet_rng( - fout.variables[dim+"_r_wet"][bin], - fout.variables[dim+"_r_wet"][bin] + fout.variables[dim+"_dr_wet"][bin] - ) + micro.diag_wet_rng( + fout.variables[dim+"_r_wet"][bin], + fout.variables[dim+"_r_wet"][bin] + fout.variables[dim+"_dr_wet"][bin] + ) elif dct["drwt"] == 'dry': - micro.diag_dry_rng( - fout.variables[dim+"_r_dry"][bin], - fout.variables[dim+"_r_dry"][bin] + fout.variables[dim+"_dr_dry"][bin] - ) + micro.diag_dry_rng( + fout.variables[dim+"_r_dry"][bin], + fout.variables[dim+"_r_dry"][bin] + fout.variables[dim+"_dr_dry"][bin] + ) + else: raise Exception("drwt should be wet or dry") + for vm in dct["moms"]: if type(vm) == int: # calculating moments @@ -175,6 +178,7 @@ def _output_bins(fout, t, micro, opts, spectra): micro.diag_dry_mom(vm) else: raise Exception("drwt should be wet or dry") fout.variables[dim+'_m'+str(vm)][int(t), int(bin)] = np.frombuffer(micro.outbuf()) + else: # calculate chemistry micro.diag_chem(_Chem_a_id[vm]) @@ -214,9 +218,14 @@ def _output_init(micro, opts, spectra): fout.createVariable(name+'_'+vm, 'd', ('t',name)) fout.variables[name+'_'+vm].unit = 'kg of chem species dissolved in cloud droplets (kg of dry air)^-1' else: - assert(type(vm)==int) - fout.createVariable(name+'_m'+str(vm), 'd', ('t',name)) - fout.variables[name+'_m'+str(vm)].unit = 'm^'+str(vm)+' (kg of dry air)^-1' + + assert(type(vm)==int) + fout.createVariable(name+'_m'+str(vm), 'd', ('t',name)) + fout.variables[name+'_m'+str(vm)].unit = 'm^'+str(vm)+' (kg of dry air)^-1' + fout.createVariable('number_of_rc_m0', 'd', ('t')) + fout.createVariable('number_of_rc_m1', 'd', ('t')) + fout.variables['number_of_rc_m0'].unit = 'm^'+str(vm)+' (kg of dry air)^-1' + fout.variables['number_of_rc_m1'].unit = 'm^'+str(vm)+' (kg of dry air)^-1' units = {"z" : "m", "t" : "s", "r_v" : "kg/kg", "th_d" : "K", "rhod" : "kg/m3", "p" : "Pa", "T" : "K", "RH" : "1" @@ -242,6 +251,12 @@ def _save_attrs(fout, dictnr): setattr(fout, var, val) def _output(fout, opts, micro, state, rec, spectra): + micro.diag_rw_ge_rc() + micro.diag_wet_mom(0) + fout.variables['number_of_rc_m0'][int(rec)] = np.frombuffer(micro.outbuf()) + micro.diag_rw_ge_rc() + micro.diag_wet_mom(1) + fout.variables['number_of_rc_m1'][int(rec)] = np.frombuffer(micro.outbuf()) _output_bins(fout, rec, micro, opts, spectra) _output_save(fout, state, rec) @@ -336,7 +351,8 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300., args, _, _, _ = inspect.getargvalues(inspect.currentframe()) opts = dict() for k in args: - opts[k] = locals()[k] + opts[k] = locals()[k] + # parsing json specification of output spectra spectra = json.loads(opts["out_bin"]) @@ -350,7 +366,7 @@ def parcel(dt=.1, z_max=200., w=1., T_0=300., p_0=101300., r_0 = .022 # water coontent specified with RH if ((opts["r_0"] < 0) and (opts["RH_0"] >= 0)): - r_0 = common.eps * opts["RH_0"] * common.p_vs(T_0) / (p_0 - opts["RH_0"] * common.p_vs(T_0)) + r_0 = common.eps * opts["RH_0"] * common.p_vs(T_0) / (p_0 - opts["RH_0"] * common.p_vs(T_0)) # sanity checks for arguments _arguments_checking(opts, spectra, aerosol) @@ -490,8 +506,8 @@ def _arguments_checking(opts, spectra, aerosol): if gstdev <= 0: raise Exception("standard deviation should be > 0 for aerosol[" + name + "]") # necessary? - if gstdev == 1.: - raise Exception("standard deviation should be != 1 to avoid monodisperse distribution for aerosol[" + name + "]") + # if gstdev == 1.: + # raise Exception("standard deviation should be != 1 to avoid monodisperse distribution for aerosol[" + name + "]") for name, dct in spectra.items(): # TODO: check if name is valid netCDF identifier diff --git a/plots/one_simulat/init_spectrum_plot.py b/plots/one_simulat/init_spectrum_plot.py index 06ba826..97dae65 100644 --- a/plots/one_simulat/init_spectrum_plot.py +++ b/plots/one_simulat/init_spectrum_plot.py @@ -27,9 +27,12 @@ def plot_init_spectrum(data, outfolder): import Gnuplot # size distribution parameters from Kreidenweis 2003 - n_tot = 566e6 - mean_r = 0.04e-6 - gstdev = 2 + n_tot = 125e6#, 15e6]#566e6 + mean_r = 0.011e-6#, 0.14e-6]#0.04e-6 + gstdev = 1.2#, 1.75]#2 + n_tot2 = 65e6#566e6 + mean_r2 = 0.06e-6#0.04e-6 + gstdev2 = 1.7#2 # from ncdf file attributes read out_bin parameters as a dictionary ... out_bin = eval(getattr(data, "out_bin")) @@ -49,9 +52,10 @@ def plot_init_spectrum(data, outfolder): # variables for plotting theoretical solution radii = np.logspace(-3, 1, 100) * 1e-6 theor = np.empty(radii.shape) + theor2 = np.empty(radii.shape) for it in range(radii.shape[0]): theor[it] = fn.log10_size_of_lnr(n_tot, mean_r, math.log(radii[it], 10), gstdev) - + theor2[it] = fn.log10_size_of_lnr(n_tot2, mean_r2, math.log(radii[it], 10), gstdev2) g = Gnuplot.Gnuplot() g('set term svg dynamic enhanced') g('reset') @@ -64,9 +68,10 @@ def plot_init_spectrum(data, outfolder): g('set yrange [0:800]') theory_r = Gnuplot.PlotItems.Data(radii * 2 * 1e6, theor * 1e-6, with_="lines", title="theory") + theory_r2 = Gnuplot.PlotItems.Data(radii * 2 * 1e6, theor2 * 1e-6, with_="lines", title="theory2") plot = Gnuplot.PlotItems.Data(rd * 2 * 1e6, model * 1e-6, with_="steps", title="model" ) - g.plot(theory_r, plot) + g.plot(theory_r, theory_r2, plot) def main(): # copy options from chem_conditions ... @@ -83,7 +88,9 @@ def main(): opts_dict['sd_conc'] = 1024 * 44 opts_dict['outfreq'] = 1 - opts_dict['out_bin'] = '{"drad": {"rght": 1e-6, "left": 1e-10, "drwt": "dry", "lnli": "log", "nbin": 26, "moms": [0,3]}}' + opts_dict['out_bin'] = '{"drad": {"rght": 2.5e-06, "left": 5e-09, "drwt": "dry", "lnli": "log", "nbin": 26, "moms": [0,1,3]}}' + + #{"rght": 2.5e-05, "moms": [0,1,3], "drwt": "wet", "nbin": 26, "lnli": "log", "left": 5e-07}} # run parcel parcel(**opts_dict) diff --git a/plots/one_simulat/profiles_plot.py b/plots/one_simulat/profiles_plot.py index 0f650e1..137bdaa 100644 --- a/plots/one_simulat/profiles_plot.py +++ b/plots/one_simulat/profiles_plot.py @@ -22,10 +22,10 @@ def plot_profiles(fnc, output_folder="../outputs"): plots[0].set_xlabel('p [hPa]') - plots[1].ticklabel_format(useOffset=False) + plots[1].ticklabel_format(useOffset=False) plots[1].set_xlabel('th_d [K]') plots[2].set_xlabel('T [K]') - plots[3].set_xlabel('kappa(rho_d :)) [kg/m3]') + plots[3].set_xlabel('kappa(rho_d :)) [kg/m3]') plots[4].set_xlabel('rv [g/kg]') plots[5].set_xlabel('RH') @@ -39,18 +39,18 @@ def plot_profiles(fnc, output_folder="../outputs"): plots[3].plot(fnc.variables["rhod"][:] , z) plots[4].plot(fnc.variables["r_v"][:] * 1000 , z) plots[5].plot( - fnc.variables["RH"][:] , z, + fnc.variables["RH"][:] , z, [fnc.variables["RH"][:].max()] * z.shape[0], z ) - + if not os.path.exists(output_folder): subprocess.call(["mkdir", output_folder]) plt.savefig(os.path.join(output_folder, "plot_profiles_onesim.svg")) def main(dt=1): - # running parcel model for different ways to solve for pressure ... + # running parcel model for different ways to solve for pressure ... outfile = "onesim_plot.nc" - parcel(dt=dt, outfreq = 10, outfile=outfile) + parcel(dt=dt, outfreq = 10, w=5, outfile=outfile) fnc = netcdf.netcdf_file(outfile) plot_profiles(fnc) fnc.close() diff --git a/plots/one_simulat/spectrum_plot.py b/plots/one_simulat/spectrum_plot.py index bdcaeaf..44ed05f 100644 --- a/plots/one_simulat/spectrum_plot.py +++ b/plots/one_simulat/spectrum_plot.py @@ -8,10 +8,13 @@ import numpy as np import pytest import subprocess - +from libcloudphxx import common +import ast +import math from parcel import parcel +import functions as fn -def plot_spectrum(data, outfolder): +def plot_spectrum(data, data2, outfolder): import Gnuplot g = Gnuplot.Gnuplot()# persist=1) @@ -23,24 +26,38 @@ def plot_spectrum(data, outfolder): rw = data.variables["wradii_r_wet"][:] * 1e6 rd = data.variables["dradii_r_dry"][:] * 1e6 + rw2 = data2.variables["wradii_r_wet"][:] * 1e6 + rd2 = data2.variables["dradii_r_dry"][:] * 1e6 + #TODO - add test if it is == to dr in netcdf - drw = np.empty(rw.shape) + drw = np.empty(rw.shape) drw[0] = rw[0] - 0 drw[1:] = (rw[1:] - rw[0:-1]) * 1e6 - drd = np.empty(rd.shape) + drd = np.empty(rd.shape) drd[0] = rd[0] - 0 drd[1:] = (rd[1:] - rd[0:-1]) * 1e6 + drw2 = np.empty(rw2.shape) + drw2[0] = rw2[0] - 0 + drw2[1:] = (rw2[1:] - rw2[0:-1]) * 1e6 + + drd2 = np.empty(rd2.shape) + drd2[0] = rd2[0] - 0 + drd2[1:] = (rd2[1:] - rd2[0:-1]) * 1e6 + for t in range(data.variables['t'].shape[0]): + g('set term svg dynamic enhanced') g('reset') g('set output "' + outfolder + 'plot_spec_' + str("%03d" % t) + '.svg"') - g('set logscale xy') - g('set ylabel "[mg^{-1} μm^{-1}]"') - g('set yrange [' + str(ymin) + ':' + str(ymax) + ']') + # g('set logscale xy') + # g('set ylabel "[mg^{-1} μm^{-1}]"') + # g('set yrange [' + str(ymin) + ':' + str(ymax) + ']') + g('set yrange [1:200]') + g('set xrange [0:15]') g('set grid') - g('set nokey') + # g('set nokey') # FSSP range g('set arrow from .5,' + str(ymin) + 'to .5,' + str(ymax) + 'nohead') @@ -50,29 +67,107 @@ def plot_spectrum(data, outfolder): nw = data.variables['wradii_m0'][t,:] / drw nd = data.variables['dradii_m0'][t,:] / drd - - plot_rw = Gnuplot.PlotItems.Data(rw, nw, with_="fsteps", title="wet radius") + + plot_rw = Gnuplot.PlotItems.Data(rw, nw, with_="fsteps", title="1") plot_rd = Gnuplot.PlotItems.Data(rd, nd, with_="fsteps", title="dry radius") - g.plot(plot_rw, plot_rd) + nw2 = data2.variables['wradii_m0'][t,:] / drw2 + nd2 = data2.variables['dradii_m0'][t,:] / drd2 + + plot_rw2 = Gnuplot.PlotItems.Data(rw2, nw2, with_="fsteps", title="1/10") + plot_rd2 = Gnuplot.PlotItems.Data(rd2, nd2, with_="fsteps", title="dry radius") + + + + g.plot(plot_rw, plot_rw2 , plot_rd) + +def plot_init_spectrum(data, outfolder): + """ + Plot the initial dry diameter distribution and compare it with the analitycal solution + + """ + import Gnuplot + + # size distribution parameters from Kreidenweis 2003 + # n_tot = 90e6 + # mean_r = 0.03e-6 + # gstdev = 1.28 + # n_tot2 = 15e6 + # mean_r2 = 0.14e-6 + # gstdev2 = 1.75 + + n_tot = 125e6 + mean_r = 0.011e-6 + gstdev = 1.2 + n_tot2 = 65e6 + mean_r2 = 0.06e-6 + gstdev2 = 1.7 + + + # from ncdf file attributes read out_bin parameters as a dictionary ... + out_bin = ast.literal_eval(getattr(data, "out_bin")) + # ... and check if the spacing used in the test was logarithmic + assert out_bin["wradii"]["lnli"] == 'log', "this test should be run with logarithmic spacing of bins" + + # parcel initial condition + rd = data.variables["wradii_r_wet"][:] # left bin edges + + # for comparison, model solution needs to be divided by log(d2) - log(d2) + # since the test is run with log spacing of bins log(d2) - log(d1) = const + d_log_rd = math.log(rd[2], 10) - math.log(rd[1], 10) + + # initial size distribution from the model + model = data.variables['wradii_m0'][0,:] * data.variables["rhod"][0] / d_log_rd + + # variables for plotting theoretical solution + radii = np.logspace(-3, 1, 100) * 1e-6 + theor = np.empty(radii.shape) + theor2 = np.empty(radii.shape) + for it in range(radii.shape[0]): + theor[it] = fn.log10_size_of_lnr(n_tot, mean_r, math.log(radii[it], 10), gstdev) + theor2[it] = fn.log10_size_of_lnr(n_tot2, mean_r2, math.log(radii[it], 10), gstdev2) + g = Gnuplot.Gnuplot() + g('set term svg dynamic enhanced') + g('reset') + g('set output "' + outfolder + '/init_spectrum.svg" ') + g('set logscale x') + g('set xlabel "particle dry diameter [μm]" ') + g('set ylabel "dN/dlog_{10}(D) [cm^{-3} log_{10}(size interval)]"') + g('set grid') + g('set xrange [0.001:10]') + g('set yrange [0:800]') + + theory_r = Gnuplot.PlotItems.Data(radii * 2 * 1e6, theor * 1e-6, with_="lines", title="theory") + theory_r2 = Gnuplot.PlotItems.Data(radii * 2 * 1e6, theor2 * 1e-6, with_="lines", title="theory2") + plot = Gnuplot.PlotItems.Data(rd * 2 * 1e6, model * 1e-6, with_="steps", title="model" ) + + g.plot(theory_r, theory_r2, plot) def main(): + + RH_init = .98 + T_init = 298. + p_init = 100000. + r_init = common.eps * RH_init * common.p_vs(T_init) / (p_init - RH_init * common.p_vs(T_init)) outfile = "test_spectrum.nc" - out_bin = '{"wradii": {"rght": 1e-4, "left": 1e-9, "drwt": "wet", "lnli": "log", "nbin": 26, "moms": [0]},\ - "dradii": {"rght": 1e-6, "left": 1e-9, "drwt": "dry", "lnli": "log", "nbin": 26, "moms": [0]}}' + outfile2 = "test_spectrum2.nc" + out_bin = '{"wradii": {"rght": 1e-4, "left": 1e-9, "drwt": "wet", "lnli": "lin", "nbin": 100, "moms": [0]},\ + "dradii": {"rght": 1e-6, "left": 1e-9, "drwt": "dry", "lnli": "lin", "nbin": 100, "moms": [0]}}' # run parcel run! - parcel(dt = .5, sd_conc = 1024, outfreq = 40, outfile = outfile, out_bin = out_bin) + parcel(dt = 1, T_0 = T_init, p_0 = p_init, RH_0 = .98, sstp_cond =1, z_max = 200, w=2, sd_conc = 10000, outfreq = 1, aerosol = '{"ammonium_sulfate": {"kappa": 0.61, "mean_r": [0.011e-6, 0.06e-6], "gstdev": [1.2, 1.7], "n_tot": [125e6, 65e6]}}', outfile = outfile, out_bin = out_bin) + parcel(dt = 1, T_0 = T_init, p_0 = p_init, RH_0 = .98, sstp_cond =20, z_max = 200, w=2, sd_conc = 10000, outfreq = 1, aerosol = '{"ammonium_sulfate": {"kappa": 0.61, "mean_r": [0.011e-6, 0.06e-6], "gstdev": [1.2, 1.7], "n_tot": [125e6, 65e6]}}', outfile = outfile2, out_bin = out_bin) data = netcdf.netcdf_file(outfile, "r") - - # plotting - plot_spectrum(data, outfolder="../outputs/") + data2 = netcdf.netcdf_file(outfile2, "r") + # plotting + plot_spectrum(data, data2, outfolder="../outputs/") + # doing plotting + plot_init_spectrum(data, outfolder="../outputs/") # cleanup subprocess.call(["rm", outfile]) if __name__ == '__main__': main() -