Solving a gray scott-system
This example is just an extension of the 2D Gray scott example. Here we show how to solve a non-linear reaction diffusion system in 3D
In figure is the final solution of the problem
More or less this example is the adaptation of the previous example to 3D with the improvement of using stencil iterator.
Stencil iterator
Stencil iterator require that you define a stencil,
{0,0,-1},
{0,0,1},
{0,-1,0},
{0,1,0},
{-1,0,0},
{1,0,0}};
once is defined it is possible get and use a stencil iterator
auto it = Old.getDomainIteratorStencil(star_stencil_3D);
while (it.isNext())
{
auto Cp = it.getStencil<0>();
auto mx = it.getStencil<1>();
auto px = it.getStencil<2>();
auto my = it.getStencil<3>();
auto py = it.getStencil<4>();
auto mz = it.getStencil<5>();
auto pz = it.getStencil<6>();
New.get<U>(Cp) = Old.get<U>(Cp) + uFactor * (
Old.get<U>(mz) +
Old.get<U>(pz) +
Old.get<U>(my) +
Old.get<U>(py) +
Old.get<U>(mx) +
Old.get<U>(px) -
6.0*Old.get<U>(Cp)) +
- deltaT * Old.get<U>(Cp) * Old.get<V>(Cp) * Old.get<V>(Cp) +
- deltaT * F * (Old.get<U>(Cp) - 1.0);
New.get<V>(Cp) = Old.get<V>(Cp) + vFactor * (
Old.get<V>(mz) +
Old.get<V>(pz) +
Old.get<V>(my) +
Old.get<V>(py) +
Old.get<V>(mx) +
Old.get<V>(px) -
6*Old.get<V>(Cp)) +
deltaT * Old.get<U>(Cp) * Old.get<V>(Cp) * Old.get<V>(Cp) +
- deltaT * (F+K) * Old.get<V>(Cp);
++it;
}
The rest of the example remain the same with the exception that the code has been extended in 3D.
- See Also
- Solve equation
Finalize
Deinitialize the library
Full code
#include "Grid/grid_dist_id.hpp"
#include "data_type/aggregate.hpp"
#include "timer.hpp"
constexpr int U = 0;
constexpr int V = 1;
constexpr int x = 0;
constexpr int y = 1;
constexpr int z = 2;
{
auto it = Old.getDomainIterator();
while (it.isNext())
{
Old.template get<U>(
key) = 1.0;
Old.template get<V>(
key) = 0.0;
New.template get<U>(
key) = 0.0;
New.template get<V>(
key) = 0.0;
++it;
}
long int x_start = Old.size(0)*1.55f/domain.
getHigh(0);
long int y_start = Old.size(1)*1.55f/domain.
getHigh(1);
long int z_start = Old.size(1)*1.55f/domain.
getHigh(2);
long int x_stop = Old.size(0)*1.85f/domain.
getHigh(0);
long int y_stop = Old.size(1)*1.85f/domain.
getHigh(1);
long int z_stop = Old.size(1)*1.85f/domain.
getHigh(2);
auto it_init = Old.getSubDomainIterator(start,stop);
while (it_init.isNext())
{
Old.template get<U>(
key) = 0.5 + (((
double)std::rand())/RAND_MAX -0.5)/10.0;
Old.template get<V>(
key) = 0.25 + (((
double)std::rand())/RAND_MAX -0.5)/20.0;
++it_init;
}
}
int main(int argc, char* argv[])
{
openfpm_init(&argc,&argv);
size_t sz[3] = {128,128,128};
double deltaT = 1;
double du = 2*1e-5;
double dv = 1*1e-5;
size_t timeSteps = 5000;
double K = 0.053;
double F = 0.014;
double spacing[3] = {Old.
spacing(0),Old.spacing(1),Old.spacing(2)};
init(Old,New,domain);
size_t count = 0;
Old.template ghost_get<U,V>();
double uFactor = deltaT * du/(spacing[x]*spacing[x]);
double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
{0,0,-1},
{0,0,1},
{0,-1,0},
{0,1,0},
{-1,0,0},
{1,0,0}};
for (size_t i = 0; i < timeSteps; ++i)
{
if (i % 300 == 0)
{std::cout << "STEP: " << i << std::endl;}
auto it = Old.getDomainIteratorStencil(star_stencil_3D);
while (it.isNext())
{
auto Cp = it.getStencil<0>();
auto mx = it.getStencil<1>();
auto px = it.getStencil<2>();
auto my = it.getStencil<3>();
auto py = it.getStencil<4>();
auto mz = it.getStencil<5>();
auto pz = it.getStencil<6>();
New.
get<U>(Cp) = Old.get<U>(Cp) + uFactor * (
Old.get<U>(mz) +
Old.get<U>(pz) +
Old.get<U>(my) +
Old.get<U>(py) +
Old.get<U>(mx) +
Old.get<U>(px) -
6.0*Old.get<U>(Cp)) +
- deltaT * Old.get<U>(Cp) * Old.get<V>(Cp) * Old.get<V>(Cp) +
- deltaT * F * (Old.get<U>(Cp) - 1.0);
New.get<V>(Cp) = Old.get<V>(Cp) + vFactor * (
Old.get<V>(mz) +
Old.get<V>(pz) +
Old.get<V>(my) +
Old.get<V>(py) +
Old.get<V>(mx) +
Old.get<V>(px) -
6*Old.get<V>(Cp)) +
deltaT * Old.get<U>(Cp) * Old.get<V>(Cp) * Old.get<V>(Cp) +
- deltaT * (F+K) * Old.get<V>(Cp);
++it;
}
Old.copy(New);
Old.ghost_get<U,V>();
if (i % 500 == 0)
{
Old.save("output_" + std::to_string(count));
count++;
}
}
std::cout <<
"Total simulation: " << tot_sim.
getwct() << std::endl;
openfpm_finalize();
}