Skip to content
Snippets Groups Projects
Commit c509de99 authored by Niklas Schnierstein's avatar Niklas Schnierstein
Browse files

minimal ou test compiles and executes

parent 84f759d5
No related branches found
No related tags found
No related merge requests found
......@@ -24,12 +24,13 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
strncpy(this->name, NAME, 256);
this->name[255] = '\0';
this->iteration = 0;
std::cerr << "3" << std::endl;
this->ou_field = new field<rnumber,be,THREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
std::cerr << "4" << std::endl;
*this->ou_field = 0.0;
this->ou_field->dft(); //have it complex all zero TODO:Ask Chichi
std::cerr << "5" << std::endl;
this->B = new field<rnumber,be,THREExTHREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
......@@ -40,7 +41,9 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
this->ou_kmax_squ = pow(ou_kmax,2);
this->ou_energy_amplitude = ou_energy_amplitude;
assert(this->kk->kM2 >= this->ou_kmax_squ); //I can avoid one if statement later in the step method
std::cerr << "6" << std::endl;
this->initialize_B();
std::cerr << "7" << std::endl;
}
......@@ -93,7 +96,7 @@ template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
{
TIMEZONE("ornstein_uhlenbeck_process::initialize_B");
*this->B = 0.0;
this->kk->CLOOP(
[&](ptrdiff_t cindex,
ptrdiff_t xindex,
......@@ -102,43 +105,49 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
double ksqu = pow(this->kk->kx[xindex],2)+
pow(this->kk->ky[yindex],2)+
pow(this->kk->kz[zindex],2);
double kabs = sqrt(ksqu);
double sigma =
this->ou_energy_amplitude *
sqrt(4*this->gamma(kabs)*pow(kabs,-5./3.)
/this->kk->nshell[(int)(kabs/this->kk->dk)]);
for(int i=0;i<3;i++) {
for(int j=0;j<3;j++) {
if (i+j == 0)
this->B->cval(cindex,i,j,0) =
sigma/2. * (1-pow(this->kk->kx[xindex],2)/ksqu);
if (ksqu <= this->ou_kmax_squ && ksqu >= this->ou_kmin_squ)
{
double kabs = sqrt(ksqu);
if (ksqu == 0) {ksqu = 1; kabs = 1;}
if (this->kk->nshell[(int)(kabs/this->kk->dk)] == 0) std::cerr << kabs << std::endl; //DELETE
double sigma =
this->ou_energy_amplitude *
sqrt(4*this->gamma(ksqu)*pow(kabs,-5./3.)
/this->kk->nshell[(int)(kabs/this->kk->dk)]);
for(int i=0;i<3;i++) {
for(int j=0;j<3;j++) {
if (i+j == 0)
this->B->cval(cindex,i,j,0) =
sigma/2. * (1-pow(this->kk->kx[xindex],2)/ksqu);
if (i+j == 4)
this->B->cval(cindex,i,j,0) =
sigma/2. * (1-pow(this->kk->kz[zindex],2)/ksqu);
if (i+j == 4)
this->B->cval(cindex,i,j,0) =
sigma/2. * (1-pow(this->kk->kz[zindex],2)/ksqu);
if (i+j == 1)
this->B->cval(cindex,i,j,0) =
sigma/2. * (0-this->kk->kx[xindex]*this->kk->ky[yindex]/ksqu);
if (i+j == 1)
this->B->cval(cindex,i,j,0) =
sigma/2. * (0-this->kk->kx[xindex]*this->kk->ky[yindex]/ksqu);
if (i+j == 3)
this->B->cval(cindex,i,j,0) =
sigma/2. * (0-this->kk->kz[zindex]*this->kk->ky[yindex]/ksqu);
if (i+j == 3)
this->B->cval(cindex,i,j,0) =
sigma/2. * (0-this->kk->kz[zindex]*this->kk->ky[yindex]/ksqu);
if (i+j == 2) {
if (i+j == 2) {
if(i==j)
this->B->cval(cindex,i,j,0) =
sigma/2. * (1-pow(this->kk->ky[yindex],2)/ksqu);
if(i==j)
this->B->cval(cindex,i,j,0) =
sigma/2. * (1-pow(this->kk->ky[yindex],2)/ksqu);
if(i!=j)
this->B->cval(cindex,i,j,0) =
sigma/2. * (0-this->kk->kx[xindex]*this->kk->kz[zindex]/ksqu);
if(i!=j)
this->B->cval(cindex,i,j,0) =
sigma/2. * (0-this->kk->kx[xindex]*this->kk->kz[zindex]/ksqu);
}
}
}
}
});
}
......
......@@ -9,15 +9,16 @@ template <typename rnumber>
int ornstein_uhlenbeck_test<rnumber>::initialize(void)
{
TIMEZONE("ornstein_uhlenbeck_test::initialize");
std::cerr << "0" << std::endl;
this->read_parameters();
std::cerr << "1" << std::endl;
this->ou = new ornstein_uhlenbeck_process<rnumber,FFTW>(
"test",
this->nx,this->ny,this->nz,
this->ou_kmin, this->ou_kmax, this->ou_energy_amplitude,
this->dkx,this->dky,this->dkz,
FFTW_ESTIMATE);
std::cerr << "2" << std::endl;
if (this->myrank == 0)
{
hid_t stat_file = H5Fopen(
......@@ -47,9 +48,14 @@ int ornstein_uhlenbeck_test<rnumber>::read_parameters()
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
H5P_DEFAULT);
this->ou_kmin = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_kmin");
this->ou_kmax = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_kmax");
this->ou_energy_amplitude = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_energy_amplitude");
// this->ou_kmin = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_kmin");
// this->ou_kmax = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_kmax");
// this->ou_energy_amplitude = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_energy_amplitude");
//TODO
this->ou_kmin = 0;
this->ou_kmax = 10;
this->ou_energy_amplitude = 1;
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
......@@ -63,6 +69,7 @@ int ornstein_uhlenbeck_test<rnumber>::do_work(void)
{
this->ou->step(1e-3);
}
std::cerr << "hello" << std::endl;
return EXIT_SUCCESS;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment