diff --git a/examples/field_convergence/EDNS.py b/examples/field_convergence/EDNS.py
index 8f391ad9aa3351b9e58a89067b6cc4d5e58f7d39..e120fc9438da541afa80bf6ff85c2038032ea8de 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 2488cfb799e2d3b6da6c96542f58b33223ce62cf..735ed4f747f4d2f25ba0418427a998bb6a26f670 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 154818658f01933b99a4c5f5a9cd4aee7c4a7eb8..b7e50d6098bf442c51e268157f68736beb7eafff 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):