diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index 37ea946e4a61151a202ced69d8a33945b3253516..5276729defc5a01d0c7415f33044f7a0bd3e2d3e 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -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);
+ }
}
}
}
+
});
}
diff --git a/cpp/full_code/ornstein_uhlenbeck_test.cpp b/cpp/full_code/ornstein_uhlenbeck_test.cpp
index d79cccb8b85a408217b3dd17f5f68f96d3512e5e..2d563bdec990d7e1da826f88eb18759ba07b5462 100644
--- a/cpp/full_code/ornstein_uhlenbeck_test.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_test.cpp
@@ -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;
}