Skip to content
Snippets Groups Projects
Commit 10b0fdbd authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

rename fields, add output

parent cab0e288
Branches
Tags
1 merge request!23WIP: Feature/use cmake
Show whitespace changes
Inline Side-by-side
bfps/cpp/full_code/symmetrize_test.cpp
+ 72
− 60
View file @ 10b0fdbd
......@@ -36,11 +36,11 @@ template <typename rnumber>
int symmetrize_test<rnumber>::do_work(void)
{
// allocate
field<rnumber, FFTW, THREE> *scal_field0 = new field<rnumber, FFTW, THREE>(
field<rnumber, FFTW, THREE> *test_field0 = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
DEFAULT_FFTW_FLAG);
field<rnumber, FFTW, THREE> *scal_field1 = new field<rnumber, FFTW, THREE>(
field<rnumber, FFTW, THREE> *test_field1 = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
DEFAULT_FFTW_FLAG);
......@@ -48,7 +48,7 @@ int symmetrize_test<rnumber>::do_work(void)
std::normal_distribution<rnumber> rdist;
rgen.seed(1);
kspace<FFTW,SMOOTH> *kk = new kspace<FFTW, SMOOTH>(
scal_field0->clayout, this->dkx, this->dky, this->dkz);
test_field0->clayout, this->dkx, this->dky, this->dkz);
if (this->myrank == 0)
{
......@@ -60,62 +60,62 @@ int symmetrize_test<rnumber>::do_work(void)
H5Fclose(stat_file);
}
// fill up scal_field0
*scal_field0 = 0.0;
scal_field0->real_space_representation = false;
// fill up test_field0
*test_field0 = 0.0;
test_field0->real_space_representation = false;
kk->CLOOP_K2(
[&](ptrdiff_t cindex,
ptrdiff_t xindex,
ptrdiff_t yindex,
ptrdiff_t zindex,
double k2){
//if (k2 < kk->kM2)
if (k2 < kk->kM2)
{
scal_field0->cval(cindex, 0, 0) = rdist(rgen);
scal_field0->cval(cindex, 0, 1) = rdist(rgen);
scal_field0->cval(cindex, 1, 0) = rdist(rgen);
scal_field0->cval(cindex, 1, 1) = rdist(rgen);
scal_field0->cval(cindex, 2, 0) = rdist(rgen);
scal_field0->cval(cindex, 2, 1) = rdist(rgen);
test_field0->cval(cindex, 0, 0) = rdist(rgen);
test_field0->cval(cindex, 0, 1) = rdist(rgen);
test_field0->cval(cindex, 1, 0) = rdist(rgen);
test_field0->cval(cindex, 1, 1) = rdist(rgen);
test_field0->cval(cindex, 2, 0) = rdist(rgen);
test_field0->cval(cindex, 2, 1) = rdist(rgen);
}
if (k2 > 0)
{
scal_field0->cval(cindex, 0, 0) /= sqrt(k2);
scal_field0->cval(cindex, 0, 1) /= sqrt(k2);
scal_field0->cval(cindex, 1, 0) /= sqrt(k2);
scal_field0->cval(cindex, 1, 1) /= sqrt(k2);
scal_field0->cval(cindex, 2, 0) /= sqrt(k2);
scal_field0->cval(cindex, 2, 1) /= sqrt(k2);
test_field0->cval(cindex, 0, 0) /= sqrt(k2);
test_field0->cval(cindex, 0, 1) /= sqrt(k2);
test_field0->cval(cindex, 1, 0) /= sqrt(k2);
test_field0->cval(cindex, 1, 1) /= sqrt(k2);
test_field0->cval(cindex, 2, 0) /= sqrt(k2);
test_field0->cval(cindex, 2, 1) /= sqrt(k2);
}
else
{
scal_field0->cval(cindex, 0, 0) = 0;
scal_field0->cval(cindex, 0, 1) = 0;
scal_field0->cval(cindex, 1, 0) = 0;
scal_field0->cval(cindex, 1, 1) = 0;
scal_field0->cval(cindex, 2, 0) = 0;
scal_field0->cval(cindex, 2, 1) = 0;
test_field0->cval(cindex, 0, 0) = 0;
test_field0->cval(cindex, 0, 1) = 0;
test_field0->cval(cindex, 1, 0) = 0;
test_field0->cval(cindex, 1, 1) = 0;
test_field0->cval(cindex, 2, 0) = 0;
test_field0->cval(cindex, 2, 1) = 0;
}
});
// dealias (?!)
//kk->template low_pass<rnumber, THREE>(scal_field0->get_cdata(), kk->kM);
kk->template dealias<rnumber, THREE>(scal_field0->get_cdata());
//kk->template low_pass<rnumber, THREE>(test_field0->get_cdata(), kk->kM);
kk->template dealias<rnumber, THREE>(test_field0->get_cdata());
// make the field divergence free
kk->template force_divfree<rnumber>(scal_field0->get_cdata());
// apply symmetrize to scal_field0
scal_field0->symmetrize();
kk->template force_divfree<rnumber>(test_field0->get_cdata());
// apply symmetrize to test_field0
//test_field0->symmetrize();
// make copy in scal_field1
// this MUST be made after symmetrizing scal_field0
// (alternatively, we may symmetrize scal_field1 as well before the ift-dft cycle
scal_field1->real_space_representation = false;
*scal_field1 = scal_field0->get_cdata();
// make copy in test_field1
// this MUST be made after symmetrizing test_field0
// (alternatively, we may symmetrize test_field1 as well before the ift-dft cycle
test_field1->real_space_representation = false;
*test_field1 = test_field0->get_cdata();
// go back and forth with scal_field1, to enforce symmetry
scal_field1->ift();
scal_field1->dft();
scal_field1->normalize();
// go back and forth with test_field1, to enforce symmetry
test_field1->ift();
test_field1->dft();
test_field1->normalize();
// now compare the two fields
double max_diff = 0;
......@@ -129,33 +129,33 @@ int symmetrize_test<rnumber>::do_work(void)
ptrdiff_t yindex,
ptrdiff_t zindex,
double k2){
double diff_re0 = scal_field0->cval(cindex, 0, 0) - scal_field1->cval(cindex, 0, 0);
double diff_re1 = scal_field0->cval(cindex, 1, 0) - scal_field1->cval(cindex, 1, 0);
double diff_re2 = scal_field0->cval(cindex, 2, 0) - scal_field1->cval(cindex, 2, 0);
double diff_im0 = scal_field0->cval(cindex, 0, 1) - scal_field1->cval(cindex, 0, 1);
double diff_im1 = scal_field0->cval(cindex, 1, 1) - scal_field1->cval(cindex, 1, 1);
double diff_im2 = scal_field0->cval(cindex, 2, 1) - scal_field1->cval(cindex, 2, 1);
double diff_re0 = test_field0->cval(cindex, 0, 0) - test_field1->cval(cindex, 0, 0);
double diff_re1 = test_field0->cval(cindex, 1, 0) - test_field1->cval(cindex, 1, 0);
double diff_re2 = test_field0->cval(cindex, 2, 0) - test_field1->cval(cindex, 2, 0);
double diff_im0 = test_field0->cval(cindex, 0, 1) - test_field1->cval(cindex, 0, 1);
double diff_im1 = test_field0->cval(cindex, 1, 1) - test_field1->cval(cindex, 1, 1);
double diff_im2 = test_field0->cval(cindex, 2, 1) - test_field1->cval(cindex, 2, 1);
double diff = sqrt(diff_re0*diff_re0 + diff_re1*diff_re1 + diff_re2*diff_re2 +
diff_im0*diff_im0 + diff_im1*diff_im1 + diff_im2*diff_im2);
double amplitude0 = (scal_field0->cval(cindex, 0, 0)*scal_field0->cval(cindex, 0, 0) +
scal_field0->cval(cindex, 1, 0)*scal_field0->cval(cindex, 1, 0) +
scal_field0->cval(cindex, 2, 0)*scal_field0->cval(cindex, 2, 0) +
scal_field0->cval(cindex, 0, 1)*scal_field0->cval(cindex, 0, 1) +
scal_field0->cval(cindex, 1, 1)*scal_field0->cval(cindex, 1, 1) +
scal_field0->cval(cindex, 2, 1)*scal_field0->cval(cindex, 2, 1));
double amplitude1 = (scal_field1->cval(cindex, 0, 0)*scal_field1->cval(cindex, 0, 0) +
scal_field1->cval(cindex, 1, 0)*scal_field1->cval(cindex, 1, 0) +
scal_field1->cval(cindex, 2, 0)*scal_field1->cval(cindex, 2, 0) +
scal_field1->cval(cindex, 0, 1)*scal_field1->cval(cindex, 0, 1) +
scal_field1->cval(cindex, 1, 1)*scal_field1->cval(cindex, 1, 1) +
scal_field1->cval(cindex, 2, 1)*scal_field1->cval(cindex, 2, 1));
double amplitude0 = (test_field0->cval(cindex, 0, 0)*test_field0->cval(cindex, 0, 0) +
test_field0->cval(cindex, 1, 0)*test_field0->cval(cindex, 1, 0) +
test_field0->cval(cindex, 2, 0)*test_field0->cval(cindex, 2, 0) +
test_field0->cval(cindex, 0, 1)*test_field0->cval(cindex, 0, 1) +
test_field0->cval(cindex, 1, 1)*test_field0->cval(cindex, 1, 1) +
test_field0->cval(cindex, 2, 1)*test_field0->cval(cindex, 2, 1));
double amplitude1 = (test_field1->cval(cindex, 0, 0)*test_field1->cval(cindex, 0, 0) +
test_field1->cval(cindex, 1, 0)*test_field1->cval(cindex, 1, 0) +
test_field1->cval(cindex, 2, 0)*test_field1->cval(cindex, 2, 0) +
test_field1->cval(cindex, 0, 1)*test_field1->cval(cindex, 0, 1) +
test_field1->cval(cindex, 1, 1)*test_field1->cval(cindex, 1, 1) +
test_field1->cval(cindex, 2, 1)*test_field1->cval(cindex, 2, 1));
double amplitude = sqrt((amplitude0 + amplitude1)/2);
if (amplitude > 0)
if (diff/amplitude > max_diff)
{
max_diff = diff / amplitude;
ix = xindex;
iy = yindex + scal_field0->clayout->starts[0];
iy = yindex + test_field0->clayout->starts[0];
iz = zindex;
k_at_max_diff = sqrt(k2);
a0 = sqrt(amplitude0);
......@@ -165,10 +165,22 @@ int symmetrize_test<rnumber>::do_work(void)
DEBUG_MSG("found maximum relative difference %g at ix = %ld, iy = %ld, iz = %ld, wavenumber = %g, amplitudes %g %g\n",
max_diff, ix, iy, iz, k_at_max_diff, a0, a1);
test_field1->io(
this->simname + "_fields.h5",
"field1",
0,
false);
test_field1->ift();
test_field1->io(
this->simname + "_fields.h5",
"field1",
0,
false);
// deallocate
delete kk;
delete scal_field1;
delete scal_field0;
delete test_field1;
delete test_field0;
return EXIT_SUCCESS;
}
......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment