diff --git a/bfps/DNS.py b/bfps/DNS.py
index 48d5c455a70f2adbf374979925e46790b031d519..32b559e12beef1ca8bd21632d3aa1ba31a4ad9ac 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -151,6 +151,7 @@ class DNS(_code):
self.parameters['nu'] = float(0.1)
self.parameters['fmode'] = int(1)
self.parameters['famplitude'] = float(0.5)
+ self.parameters['friction_coefficient'] = float(0.5)
self.parameters['energy'] = float(0.5)
self.parameters['injection_rate'] = float(0.4)
self.parameters['fk0'] = float(2.0)
diff --git a/bfps/cpp/full_code/NSVE.cpp b/bfps/cpp/full_code/NSVE.cpp
index ba8a3ed65eaf17738e363f3e8c82246f1dfa5dd7..e8bf9fd2e9786bd5df8a4c9d7dc9b37e15de85c1 100644
--- a/bfps/cpp/full_code/NSVE.cpp
+++ b/bfps/cpp/full_code/NSVE.cpp
@@ -47,6 +47,7 @@ int NSVE<rnumber>::initialize(void)
this->fs->nu = nu;
this->fs->fmode = fmode;
this->fs->famplitude = famplitude;
+ this->fs->friction_coefficient = friction_coefficient;
this->fs->energy = energy;
this->fs->injection_rate = injection_rate;
this->fs->fk0 = fk0;
diff --git a/bfps/cpp/full_code/NSVE.hpp b/bfps/cpp/full_code/NSVE.hpp
index e3f6b2765875084ed152e914cd285e331502d673..062627fd1cc9513bbd29a14199e05d3a084c0851 100644
--- a/bfps/cpp/full_code/NSVE.hpp
+++ b/bfps/cpp/full_code/NSVE.hpp
@@ -42,6 +42,7 @@ class NSVE: public direct_numerical_simulation
/* parameters that are read in read_parameters */
double dt;
double famplitude;
+ double friction_coefficient;
double fk0;
double fk1;
double energy;
diff --git a/bfps/cpp/full_code/postprocess.cpp b/bfps/cpp/full_code/postprocess.cpp
index edb5929f72c5197c123f8f4e20d426ca1ad9eb6f..c48bcdb85ae496cf9ea606a20b9434c3c5abe99e 100644
--- a/bfps/cpp/full_code/postprocess.cpp
+++ b/bfps/cpp/full_code/postprocess.cpp
@@ -57,6 +57,9 @@ int postprocess::read_parameters()
dset = H5Dopen(parameter_file, "/parameters/famplitude", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->famplitude);
H5Dclose(dset);
+ dset = H5Dopen(parameter_file, "/parameters/friction_coefficient", H5P_DEFAULT);
+ H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->friction_coefficient);
+ H5Dclose(dset);
dset = H5Dopen(parameter_file, "/parameters/fk0", H5P_DEFAULT);
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->fk0);
H5Dclose(dset);
diff --git a/bfps/cpp/full_code/postprocess.hpp b/bfps/cpp/full_code/postprocess.hpp
index c80fc3f2dfdc35691d9e69442fa3ad7b6e592891..660e561a1ddba6b355ee5ecf49c5cea1b7153eff 100644
--- a/bfps/cpp/full_code/postprocess.hpp
+++ b/bfps/cpp/full_code/postprocess.hpp
@@ -43,6 +43,7 @@ class postprocess: public code_base
/* parameters that are read in read_parameters */
double dt;
double famplitude;
+ double friction_coefficient;
double fk0;
double fk1;
int fmode;
diff --git a/bfps/cpp/vorticity_equation.cpp b/bfps/cpp/vorticity_equation.cpp
index 5c1cc056282ba9a1b490e6f8a9ae423bc880faab..47fb0c68907e441e9b49df7044755262b7f296be 100644
--- a/bfps/cpp/vorticity_equation.cpp
+++ b/bfps/cpp/vorticity_equation.cpp
@@ -152,6 +152,7 @@ vorticity_equation<rnumber, be>::vorticity_equation(
this->nu = 0.1;
this->fmode = 1;
this->famplitude = 1.0;
+ this->friction_coefficient = 1.0;
this->fk0 = 2.0;
this->fk1 = 4.0;
}
@@ -225,6 +226,53 @@ void vorticity_equation<rnumber, be>::compute_velocity(field<rnumber, be, THREE>
this->u->symmetrize();
}
+template <class rnumber,
+ field_backend be>
+void vorticity_equation<rnumber, be>::add_Kolmogorov_forcing(
+ field<rnumber, be, THREE> *dst,
+ int fmode,
+ double famplitude)
+{
+ TIMEZONE("vorticity_equation::add_Kolmogorov_forcing");
+ ptrdiff_t cindex;
+ if (dst->clayout->myrank == dst->clayout->rank[0][fmode])
+ {
+ cindex = dst->get_cindex(0, (fmode - dst->clayout->starts[0]), 0);
+ dst->cval(cindex,2, 0) -= famplitude/2;
+ }
+ if (dst->clayout->myrank == dst->clayout->rank[0][dst->clayout->sizes[0] - fmode])
+ {
+ cindex = dst->get_cindex(0, (dst->clayout->sizes[0] - fmode - dst->clayout->starts[0]), 0);
+ dst->cval(cindex, 2, 0) -= famplitude/2;
+ }
+}
+
+template <class rnumber,
+ field_backend be>
+void vorticity_equation<rnumber, be>::add_field_band(
+ field<rnumber, be, THREE> *dst,
+ field<rnumber, be, THREE> *src,
+ double k0, double k1,
+ double prefactor)
+{
+ TIMEZONE("vorticity_equation::add_field_band");
+ this->kk->CLOOP(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex){
+ double knorm = sqrt(this->kk->kx[xindex]*this->kk->kx[xindex] +
+ this->kk->ky[yindex]*this->kk->ky[yindex] +
+ this->kk->kz[zindex]*this->kk->kz[zindex]);
+ if ((k0 <= knorm) &&
+ (k1 >= knorm))
+ for (int c=0; c<3; c++)
+ for (int i=0; i<2; i++)
+ dst->cval(cindex,c,i) += prefactor*src->cval(cindex,c,i);
+ }
+ );
+}
+
template <class rnumber,
field_backend be>
void vorticity_equation<rnumber, be>::add_forcing(
@@ -236,69 +284,37 @@ void vorticity_equation<rnumber, be>::add_forcing(
return;
if (strcmp(this->forcing_type, "Kolmogorov") == 0)
{
- ptrdiff_t cindex;
- if (this->cvorticity->clayout->myrank == this->cvorticity->clayout->rank[0][this->fmode])
- {
- cindex = dst->get_cindex(0, (this->fmode - this->cvorticity->clayout->starts[0]), 0);
- dst->cval(cindex,2, 0) -= this->famplitude/2;
- }
- if (this->cvorticity->clayout->myrank == this->cvorticity->clayout->rank[0][this->cvorticity->clayout->sizes[0] - this->fmode])
- {
- cindex = dst->get_cindex(0, (this->cvorticity->clayout->sizes[0] - this->fmode - this->cvorticity->clayout->starts[0]), 0);
- dst->cval(cindex, 2, 0) -= this->famplitude/2;
- }
+ this->add_Kolmogorov_forcing(dst, this->fmode, this->famplitude);
return;
}
if (strcmp(this->forcing_type, "2Kolmogorov") == 0)
{
// 2 Kolmogorov forces
// first one wavenumber fk0, amplitude 1 - A
- DEBUG_MSG("famplitude = %g\n", this->famplitude);
- ptrdiff_t cindex;
double amplitude = 1 - this->famplitude;
int fmode = int(this->fk0 / this->kk->dky);
- if (this->cvorticity->clayout->myrank == this->cvorticity->clayout->rank[0][fmode])
- {
- cindex = dst->get_cindex(0, (fmode - this->cvorticity->clayout->starts[0]), 0);
- dst->cval(cindex,2, 0) -= amplitude/2;
- }
- if (this->cvorticity->clayout->myrank == this->cvorticity->clayout->rank[0][this->cvorticity->clayout->sizes[0] - fmode])
- {
- cindex = dst->get_cindex(0, (this->cvorticity->clayout->sizes[0] - fmode - this->cvorticity->clayout->starts[0]), 0);
- dst->cval(cindex, 2, 0) -= amplitude/2;
- }
+ this->add_Kolmogorov_forcing(dst, fmode, amplitude);
// second one wavenumber fk1, amplitude A
amplitude = this->famplitude * pow(int(this->fk1) / double(int(this->fk0)), 3);
fmode = int(this->fk1 / this->kk->dky);
- if (this->cvorticity->clayout->myrank == this->cvorticity->clayout->rank[0][fmode])
- {
- cindex = dst->get_cindex(0, (fmode - this->cvorticity->clayout->starts[0]), 0);
- dst->cval(cindex,2, 0) -= amplitude/2;
- }
- if (this->cvorticity->clayout->myrank == this->cvorticity->clayout->rank[0][this->cvorticity->clayout->sizes[0] - fmode])
- {
- cindex = dst->get_cindex(0, (this->cvorticity->clayout->sizes[0] - fmode - this->cvorticity->clayout->starts[0]), 0);
- dst->cval(cindex, 2, 0) -= amplitude/2;
- }
+ this->add_Kolmogorov_forcing(dst, fmode, amplitude);
+ return;
+ }
+ if (strcmp(this->forcing_type, "Kolmogorov_and_drag") == 0)
+ {
+ this->add_Kolmogorov_forcing(dst, this->fmode, this->famplitude);
+ this->add_field_band(
+ dst, vort_field,
+ 0, this->fmode,
+ -this->friction_coefficient);
return;
}
if (strcmp(this->forcing_type, "linear") == 0)
{
- this->kk->CLOOP(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
- double knorm = sqrt(this->kk->kx[xindex]*this->kk->kx[xindex] +
- this->kk->ky[yindex]*this->kk->ky[yindex] +
- this->kk->kz[zindex]*this->kk->kz[zindex]);
- if ((this->fk0 <= knorm) &&
- (this->fk1 >= knorm))
- for (int c=0; c<3; c++)
- for (int i=0; i<2; i++)
- dst->cval(cindex,c,i) += this->famplitude*vort_field->cval(cindex,c,i);
- }
- );
+ this->add_field_band(
+ dst, vort_field,
+ this->fk0, this->fk1,
+ this->famplitude);
return;
}
if (strcmp(this->forcing_type, "fixed_energy_injection_rate") == 0)
@@ -339,20 +355,10 @@ void vorticity_equation<rnumber, be>::add_forcing(
// now, modify amplitudes
double temp_famplitude = this->injection_rate / energy_in_shell;
- this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
- double knorm = sqrt(k2);
- if ((this->fk0 <= knorm) &&
- (this->fk1 >= knorm))
- for (int c=0; c<3; c++)
- for (int i=0; i<2; i++)
- dst->cval(cindex,c,i) += temp_famplitude*vort_field->cval(cindex, c, i);
- }
- );
+ this->add_field_band(
+ dst, vort_field,
+ this->fk0, this->fk1,
+ temp_famplitude);
return;
}
if (strcmp(this->forcing_type, "fixed_energy") == 0)
diff --git a/bfps/cpp/vorticity_equation.hpp b/bfps/cpp/vorticity_equation.hpp
index 7ce310712323b7c943d60c75a3f9f45606c1caa1..81f0cb667e5377ae34ee90ef94e62c45b7acc541 100644
--- a/bfps/cpp/vorticity_equation.hpp
+++ b/bfps/cpp/vorticity_equation.hpp
@@ -67,11 +67,12 @@ class vorticity_equation
/* physical parameters */
double nu;
- int fmode; // for Kolmogorov flow
- double famplitude; // both for Kflow and band forcing
- double fk0, fk1; // for band forcing
- double injection_rate; // for fixed energy injection rate
- double energy; // for fixed energy
+ int fmode; // for Kolmogorov flow
+ double famplitude; // both for Kflow and band forcing
+ double fk0, fk1; // for band forcing
+ double injection_rate; // for fixed energy injection rate
+ double energy; // for fixed energy
+ double friction_coefficient; // for Kolmogorov_and_drag
char forcing_type[128];
/* constructor, destructor */
@@ -101,6 +102,15 @@ class vorticity_equation
void add_forcing(field<rnumber, be, THREE> *dst,
field<rnumber, be, THREE> *src_vorticity);
+ void add_Kolmogorov_forcing(field<rnumber, be, THREE> *dst,
+ int fmode,
+ double famplitude);
+ void add_field_band(
+ field<rnumber, be, THREE> *dst,
+ field<rnumber, be, THREE> *src,
+ double k0, double k1,
+ double prefactor);
+
/** \brief Method that imposes action of forcing on new vorticity field.
*
* If the force is implicit, in the sense that kinetic energy must be