OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
Gray Scott in 3D

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,

static grid_key_dx<3> star_stencil_3D[7] = {{0,0,0},
{0,0,-1},
{0,0,1},
{0,-1,0},
{0,1,0},
{-1,0,0},
{1,0,0}};
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19

once is defined it is possible get and use a stencil iterator

auto it = Old.getDomainIteratorStencil(star_stencil_3D);
while (it.isNext())
{
// center point
auto Cp = it.getStencil<0>();
// plus,minus X,Y,Z
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>();
// update based on Eq 2
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);
// update based on Eq 2
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);
// Next point in the grid
++it;
}
[v_transform metafunction]

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

openfpm_finalize();

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;
void init(grid_dist_id<3,double,aggregate<double,double> > & Old, grid_dist_id<3,double,aggregate<double,double> > & New, Box<3,double> & domain)
{
auto it = Old.getDomainIterator();
while (it.isNext())
{
// Get the local grid key
auto key = it.get();
// Old values U and V
Old.template get<U>(key) = 1.0;
Old.template get<V>(key) = 0.0;
// Old values U and V
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);
grid_key_dx<3> start({x_start,y_start,z_start});
grid_key_dx<3> stop ({x_stop,y_stop,z_stop});
auto it_init = Old.getSubDomainIterator(start,stop);
while (it_init.isNext())
{
auto key = it_init.get();
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);
// domain
Box<3,double> domain({0.0,0.0,0.0},{2.5,2.5,2.5});
// grid size
size_t sz[3] = {128,128,128};
// Define periodicity of the grid
periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
// Ghost in grid unit
// deltaT
double deltaT = 1;
// Diffusion constant for specie U
double du = 2*1e-5;
// Diffusion constant for specie V
double dv = 1*1e-5;
// Number of timesteps
size_t timeSteps = 5000;
// K and F (Physical constant in the equation)
double K = 0.053;
double F = 0.014;
// New grid with the decomposition of the old grid
grid_dist_id<3, double, aggregate<double,double>> New(Old.getDecomposition(),sz,g);
// spacing of the grid on x and y
double spacing[3] = {Old.spacing(0),Old.spacing(1),Old.spacing(2)};
init(Old,New,domain);
// sync the ghost
size_t count = 0;
Old.template ghost_get<U,V>();
// because we assume that spacing[x] == spacing[y] we use formula 2
// and we calculate the prefactor of Eq 2
double uFactor = deltaT * du/(spacing[x]*spacing[x]);
double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
timer tot_sim;
tot_sim.start();
static grid_key_dx<3> star_stencil_3D[7] = {{0,0,0},
{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())
{
// center point
auto Cp = it.getStencil<0>();
// plus,minus X,Y,Z
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>();
// update based on Eq 2
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);
// update based on Eq 2
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);
// Next point in the grid
++it;
}
// Here we copy New into the old grid in preparation of the new step
// It would be better to alternate, but using this we can show the usage
// of the function copy. To note that copy work only on two grid of the same
// decomposition. If you want to copy also the decomposition, or force to be
// exactly the same, use Old = New
Old.copy(New);
// After copy we synchronize again the ghost part U and V
Old.ghost_get<U,V>();
// Every 500 time step we output the configuration for
// visualization
if (i % 500 == 0)
{
// Old.save("output_" + std::to_string(count));
count++;
}
}
tot_sim.stop();
std::cout << "Total simulation: " << tot_sim.getwct() << std::endl;
openfpm_finalize();
}
This class represent an N-dimensional box.
Definition Box.hpp:61
This is a distributed grid.
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition grid_key.hpp:503
Class for cpu time benchmarking.
Definition timer.hpp:28
void stop()
Stop the timer.
Definition timer.hpp:119
void start()
Start the timer.
Definition timer.hpp:90
double getwct()
Return the elapsed real time.
Definition timer.hpp:130
OutputIteratorT OffsetT ReductionOpT OuputT init
< [in] The initial value of the reduction
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
Boundary conditions.
Definition common.hpp:22