Skip to content
GitLab
Commit e5c265c3 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/get-rfields' into develop

parents 9f533ef2 7611d440
Pipeline #65679 passed with stage
in 3 minutes and 45 seconds
Hide whitespace changes
Inline Side-by-side
TurTLE/PP.py
View file @ e5c265c3
......@@ -145,6 +145,8 @@ class PP(_code):
pars['filter_type'] = 'Gauss'
pars['max_velocity_estimate'] = float(10)
pars['histogram_bins'] = int(129)
elif dns_type == 'get_rfields':
pars['TrS2_on'] = int(0)
return pars
def get_data_file_name(self):
return os.path.join(self.work_dir, self.simname + '.h5')
......@@ -806,7 +808,7 @@ class PP(_code):
ff.require_group('vorticity')
ff.require_group('vorticity/complex')
checkpoint_file_list = [self.simname + '_checkpoint_{0}.h5'.format(cp)
for cp in range(df['checkpoint'][()])]
for cp in range(df['checkpoint'][()]+1)]
for cpf_name in checkpoint_file_list:
if os.path.exists(cpf_name):
cpf = h5py.File(cpf_name, 'r')
......
cpp/full_code/field_single_to_double.cpp
View file @ e5c265c3
......@@ -74,7 +74,9 @@ template <typename rnumber>
int field_single_to_double<rnumber>::work_on_current_iteration(void)
{
TIMEZONE("field_single_to_double::work_on_current_iteration");
DEBUG_MSG("before read_vorticity at iteration %d\n", this->iteration);
this->read_current_cvorticity();
DEBUG_MSG("after read_vorticity at iteration %d\n", this->iteration);
// using CLOOP as opposed to a global std::copy because CLOOP
// is openmp parallelized.
......
cpp/full_code/get_rfields.cpp
View file @ e5c265c3
......@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
* This file is part of bfps. *
* This file is part of TurTLE. *
* *
* bfps is free software: you can redistribute it and/or modify *
* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
* bfps is distributed in the hope that it will be useful, *
* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with bfps. If not, see <http://www.gnu.org/licenses/> *
* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
......@@ -35,8 +35,27 @@ int get_rfields<rnumber>::initialize(void)
TIMEZONE("get_rfields::initialize");
this->NSVE_field_stats<rnumber>::initialize();
DEBUG_MSG("after NSVE_field_stats::initialize\n");
// allocate kspace
this->kk = new kspace<FFTW, SMOOTH>(
this->vorticity->clayout, this->dkx, this->dky, this->dkz);
// allocate field for TrS2
this->traceS2 = new field<rnumber, FFTW, ONE>(
this->nx, this->ny, this->nz,
this->comm,
this->vorticity->fftw_plan_rigor);
this->traceS2->real_space_representation = true;
// allocate field for velocity
this->vel = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
this->vorticity->fftw_plan_rigor);
this->vel->real_space_representation = false;
// allocate field for velocity gradient
this->grad_vel = new field<rnumber, FFTW, THREExTHREE>(
this->nx, this->ny, this->nz,
this->comm,
this->vorticity->fftw_plan_rigor);
hid_t parameter_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
......@@ -53,6 +72,9 @@ int get_rfields<rnumber>::initialize(void)
this->iteration_list = hdf5_tools::read_vector<int>(
parameter_file,
"/get_rfields/parameters/iteration_list");
this->TrS2_on = hdf5_tools::read_value<int>(
parameter_file,
"/get_rfields/parameters/TrS2_on");
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
......@@ -62,53 +84,66 @@ int get_rfields<rnumber>::work_on_current_iteration(void)
{
TIMEZONE("get_rfields::work_on_current_iteration");
this->read_current_cvorticity();
field<rnumber, FFTW, THREE> *vel = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
this->vorticity->fftw_plan_rigor);
vel->real_space_representation = false;
this->kk->CLOOP_K2(
[&](ptrdiff_t cindex,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex,
double k2){
if (k2 <= this->kk->kM2 && k2 > 0)
{
vel->cval(cindex,0,0) = -(this->kk->ky[yindex]*this->vorticity->cval(cindex,2,1) - this->kk->kz[zindex]*this->vorticity->cval(cindex,1,1)) / k2;
vel->cval(cindex,0,1) = (this->kk->ky[yindex]*this->vorticity->cval(cindex,2,0) - this->kk->kz[zindex]*this->vorticity->cval(cindex,1,0)) / k2;
vel->cval(cindex,1,0) = -(this->kk->kz[zindex]*this->vorticity->cval(cindex,0,1) - this->kk->kx[xindex]*this->vorticity->cval(cindex,2,1)) / k2;
vel->cval(cindex,1,1) = (this->kk->kz[zindex]*this->vorticity->cval(cindex,0,0) - this->kk->kx[xindex]*this->vorticity->cval(cindex,2,0)) / k2;
vel->cval(cindex,2,0) = -(this->kk->kx[xindex]*this->vorticity->cval(cindex,1,1) - this->kk->ky[yindex]*this->vorticity->cval(cindex,0,1)) / k2;
vel->cval(cindex,2,1) = (this->kk->kx[xindex]*this->vorticity->cval(cindex,1,0) - this->kk->ky[yindex]*this->vorticity->cval(cindex,0,0)) / k2;
}
else
std::fill_n((rnumber*)(vel->get_cdata()+3*cindex), 6, 0.0);
}
);
vel->symmetrize();
vel->ift();
invert_curl(
this->kk,
this->vorticity,
this->vel);
/// compute velocity gradient
compute_gradient<rnumber, FFTW, THREE, THREExTHREE>(
this->kk,
this->vel,
this->grad_vel);
// compute TrS2
this->grad_vel->ift();
this->traceS2->RLOOP(
[&](ptrdiff_t rindex,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex){
rnumber AxxAxx = this->grad_vel->rval(rindex, 0, 0)*this->grad_vel->rval(rindex, 0, 0);
rnumber AyyAyy = this->grad_vel->rval(rindex, 1, 1)*this->grad_vel->rval(rindex, 1, 1);
rnumber AzzAzz = this->grad_vel->rval(rindex, 2, 2)*this->grad_vel->rval(rindex, 2, 2);
rnumber Sxy = this->grad_vel->rval(rindex, 0, 1) + this->grad_vel->rval(rindex, 1, 0);
rnumber Syz = this->grad_vel->rval(rindex, 1, 2) + this->grad_vel->rval(rindex, 2, 1);
rnumber Szx = this->grad_vel->rval(rindex, 2, 0) + this->grad_vel->rval(rindex, 0, 2);
this->traceS2->rval(rindex) = (
AxxAxx + AyyAyy + AzzAzz +
(Sxy*Sxy + Syz*Syz + Szx*Szx)/2);
});
std::string fname = (
this->simname +
std::string("_checkpoint_") +
std::to_string(this->iteration / (this->niter_out*this->checkpoints_per_file)) +
std::string(".h5"));
vel->io(
/// output velocity field
this->vel->ift();
this->vel->io(
fname,
"velocity",
this->iteration,
false);
delete vel;
/// output vorticity field
this->vorticity->ift();
this->vorticity->io(
fname,
"vorticity",
this->iteration,
false);
if (this->TrS2_on)
{
this->traceS2->io(
fname,
"TrS2",
this->iteration,
false);
}
return EXIT_SUCCESS;
}
......@@ -117,6 +152,9 @@ int get_rfields<rnumber>::finalize(void)
{
TIMEZONE("get_rfields::finalize");
delete this->kk;
delete this->traceS2;
delete this->vel;
delete this->grad_vel;
this->NSVE_field_stats<rnumber>::finalize();
return EXIT_SUCCESS;
}
......
cpp/full_code/get_rfields.hpp
View file @ e5c265c3
......@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
* This file is part of bfps. *
* This file is part of TurTLE. *
* *
* bfps is free software: you can redistribute it and/or modify *
* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
* bfps is distributed in the hope that it will be useful, *
* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with bfps. If not, see <http://www.gnu.org/licenses/> *
* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
......@@ -42,7 +42,11 @@ class get_rfields: public NSVE_field_stats<rnumber>
public:
int checkpoints_per_file;
int niter_out;
int TrS2_on;
kspace<FFTW, SMOOTH> *kk;
field<rnumber, FFTW, ONE> *traceS2;
field<rnumber, FFTW, THREE> *vel;
field<rnumber, FFTW, THREExTHREE> *grad_vel;
get_rfields(
const MPI_Comm COMMUNICATOR,
......
tests/run_all_tests.sh
View file @ e5c265c3
......@@ -8,6 +8,6 @@ turtle.test_NSVEparticles
# test postprocessing
turtle PP field_single_to_double --simname dns_nsveparticles --iter0 32 --iter1 32
turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64
turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64 --TrS2_on 1
turtle PP joint_acc_vel_stats --simname dns_nsveparticles --iter0 0 --iter1 64
turtle PP resize --simname dns_nsveparticles --new_nx 96 --new_ny 96 --new_nz 96 --new_simname dns_nsveparticles_resized
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment