Commit db01eb7d authored by Chichi Lalescu's avatar Chichi Lalescu
Browse files

add new force

parent 474bdab2
Pipeline #23861 passed with stage
in 6 minutes and 23 seconds
......@@ -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)
......
......@@ -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;
......
......@@ -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;
......
......@@ -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);
......
......@@ -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;
......
......@@ -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)
......
......@@ -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
......
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