From 36333ab00393e82c1b73c50e4ea6ad1a4f64dd12 Mon Sep 17 00:00:00 2001 From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de> Date: Sat, 2 Dec 2023 06:33:10 +0100 Subject: [PATCH] adds histogram plot --- examples/field_convergence/EDNS.py | 41 ++++++++++++++++------- examples/field_convergence/NSVE_error.cpp | 33 +++++++++++------- examples/field_convergence/NSVE_error.hpp | 3 ++ 3 files changed, 53 insertions(+), 24 deletions(-) diff --git a/examples/field_convergence/EDNS.py b/examples/field_convergence/EDNS.py index 8f391ad9..e120fc94 100644 --- a/examples/field_convergence/EDNS.py +++ b/examples/field_convergence/EDNS.py @@ -83,6 +83,8 @@ class EDNS(TurTLE.DNS): def generate_default_parameters(self): TurTLE.DNS.generate_default_parameters(self) self.extra_parameters['NSVE_error'] = { + 'max_err_estimate' : 1e-5, + 'max_err_relative_estimate' : 1e-4, } return None @@ -217,24 +219,38 @@ class EDNS(TurTLE.DNS): work_dir = self.work_dir) for i in [1, 2, 4]] return member_list + def read_error_statistics( + self): + df = self.get_data_file() + for ii in [1, 2]: + self.statistics['error/moments/err{0}'.format(ii)] = df['statistics/moments/err{0}'.format(ii)][()] + self.statistics['error/moments/err{0}_relative'.format(ii)] = df['statistics/moments/err{0}_relative'.format(ii)][()] + self.statistics['error/histograms/err{0}'.format(ii)] = df['statistics/histograms/err{0}'.format(ii)][()] + self.statistics['error/histograms/err{0}_relative'.format(ii)] = df['statistics/histograms/err{0}_relative'.format(ii)][()] + df.close() + return None def plot_error_statistics( self, fname = None, err_norm = 1.0): if type(fname) == type(None): fname = self.simname + '_err_stats' - df = self.get_data_file() - f = plt.figure() - a = f.add_subplot(211) - err1 = np.sqrt(df['statistics/moments/err1'][:, 2, 3]) / err_norm - err2 = np.sqrt(df['statistics/moments/err2'][:, 2, 3]) / err_norm - a.plot(err1) - a.plot(err2) - a = f.add_subplot(212) - err1 = np.sqrt(df['statistics/moments/err1_relative'][:, 2, 3]) - err2 = np.sqrt(df['statistics/moments/err2_relative'][:, 2, 3]) - a.plot(err1) - a.plot(err2) + f = plt.figure(figsize = (8, 8)) + a = f.add_subplot(221) + a.plot(self.statistics['error/moments/err1'][:, 2, 3]**0.5) + a.plot(self.statistics['error/moments/err2'][:, 2, 3]**0.5) + a = f.add_subplot(222) + a.plot(self.statistics['error/moments/err1_relative'][:, 2, 3]**0.5) + a.plot(self.statistics['error/moments/err2_relative'][:, 2, 3]**0.5) + a = f.add_subplot(223) + a.plot(self.statistics['error/histograms/err1'][-1, :, 3]) + a.plot(self.statistics['error/histograms/err2'][-1, :, 3]) + a.set_yscale('log') + a = f.add_subplot(224) + print(self.statistics['error/histograms/err1_relative'][-1, :, 3]) + a.plot(self.statistics['error/histograms/err1_relative'][-1, :, 3]) + a.plot(self.statistics['error/histograms/err2_relative'][-1, :, 3]) + a.set_yscale('log') f.tight_layout() f.savefig(fname + '.pdf') plt.close(f) @@ -250,6 +266,7 @@ def main(): + sys.argv[1:]) bla.read_parameters() mm = bla.run_sanity_check() + bla.read_error_statistics() bla.plot_error_statistics(err_norm = mm[0].statistics['Uint']) return None diff --git a/examples/field_convergence/NSVE_error.cpp b/examples/field_convergence/NSVE_error.cpp index 2488cfb7..735ed4f7 100644 --- a/examples/field_convergence/NSVE_error.cpp +++ b/examples/field_convergence/NSVE_error.cpp @@ -121,7 +121,9 @@ int compare_real_fields( const field<rnumber, FFTW, THREE>* b, const std::string suffix, const hid_t stat_group, - const int ii) + const int ii, + const double max_err_estimate, + const double max_err_relative_estimate) { TIMEZONE("compare_real_fields"); field<rnumber, FFTW, THREE> * err; @@ -153,13 +155,13 @@ int compare_real_fields( stat_group, "err" + suffix, ii, - 1.0); + max_err_estimate); err_relative->compute_stats( kk, stat_group, "err" + suffix + "_relative", ii, - 1.0); + max_err_relative_estimate); delete err; delete err_relative; return EXIT_SUCCESS; @@ -190,7 +192,9 @@ int NSVE_error<rnumber>::do_stats() this->nsve1->tmp_vec_field, "1", stat_group, - this->iteration / niter_stat); + this->iteration / niter_stat, + this->max_err_estimate, + this->max_err_relative_estimate); compare_real_fields( this->nsve1->fs->kk, @@ -198,7 +202,9 @@ int NSVE_error<rnumber>::do_stats() this->nsve2->tmp_vec_field, "2", stat_group, - this->iteration / niter_stat); + this->iteration / niter_stat, + this->max_err_estimate, + this->max_err_relative_estimate); if (this->myrank == 0) H5Gclose(stat_group); @@ -210,13 +216,16 @@ int NSVE_error<rnumber>::read_parameters(void) { TIMEZONE("NSVE_error::read_parameters"); this->direct_numerical_simulation::read_parameters(); - //hid_t parameter_file = H5Fopen( - // (this->ensname + ".h5").c_str(), - // H5F_ACC_RDONLY, - // H5P_DEFAULT); - //// read parameters here ... - //H5Fclose(parameter_file); - //MPI_Barrier(this->comm); // "lock" on parameter file + hid_t parameter_file = H5Fopen( + (this->simname + ".h5").c_str(), + H5F_ACC_RDONLY, + H5P_DEFAULT); + this->max_err_estimate = hdf5_tools::read_value<double>( + parameter_file, "parameters/max_err_estimate"); + this->max_err_relative_estimate = hdf5_tools::read_value<double>( + parameter_file, "parameters/max_err_relative_estimate"); + H5Fclose(parameter_file); + MPI_Barrier(this->comm); // "lock" on parameter file return EXIT_SUCCESS; } diff --git a/examples/field_convergence/NSVE_error.hpp b/examples/field_convergence/NSVE_error.hpp index 15481865..b7e50d60 100644 --- a/examples/field_convergence/NSVE_error.hpp +++ b/examples/field_convergence/NSVE_error.hpp @@ -20,6 +20,9 @@ class NSVE_error: public direct_numerical_simulation NSVE<rnumber> * nsve2; NSVE<rnumber> * nsve4; + double max_err_estimate; + double max_err_relative_estimate; + NSVE_error( const MPI_Comm COMMUNICATOR, const std::string &simname): -- GitLab