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

Solving a gray scott-system

This example show the usage of periodic grid with ghost part given in grid units to solve the following system of equations

\(\frac{\partial u}{\partial t} = D_u \nabla^{2} u - uv^2 + F(1-u)\)

\(\frac{\partial v}{\partial t} = D_v \nabla^{2} v + uv^2 - (F + k)v\)

Constants and functions

First we define convenient constants

constexpr int U = 0;
constexpr int V = 1;
constexpr int x = 0;
constexpr int y = 1;

We also define an init function. This function initialize the species U and V. In the following we are going into the detail of this function

In figure is the final solution of the problem

void init(grid_dist_id<2,double,aggregate<double,double> > & Old, grid_dist_id<2,double,aggregate<double,double> > & New, Box<2,double> & domain)
{
This class represent an N-dimensional box.
Definition Box.hpp:61
This is a distributed grid.
aggregate of properties, from a list of object if create a struct that follow the OPENFPM native stru...
}

Here we initialize for the full domain. U and V itarating across the grid points. For the calculation We are using 2 grids one Old and New. We initialize Old with the initial condition concentration of the species U = 1 over all the domain and concentration of the specie V = 0 over all the domain. While the New grid is initialized to 0

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;
}

After we initialized the full grid, we create a perturbation in the domain with different values. We do in the part of space: 1.55 < x < 1.85 and 1.55 < y < 1.85. Or more precisely on the points included in this part of space.

grid_key_dx<2> start({(long int)std::floor(Old.size(0)*1.55f/domain.getHigh(0)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(1))});
grid_key_dx<2> stop ({(long int)std::ceil (Old.size(0)*1.85f/domain.getHigh(0)),(long int)std::ceil (Old.size(1)*1.85f/domain.getHigh(1))});
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)/100.0;
Old.template get<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/200.0;
++it_init;
}
grid_key_dx is the key to access any element in the grid
Definition grid_key.hpp:19
__device__ __host__ index_type get(index_type i) const
Get the i index.
Definition grid_key.hpp:503

Initialization

Initialize the library

Create

  • A 2D box that define the domain
  • an array of 2 unsigned integer that will define the size of the grid on each dimension
  • Periodicity of the grid
  • A Ghost object that will define the extension of the ghost part for each sub-domain in grid point unit

We also define numerical and physical parameters

  • Time stepping for the integration
  • Diffusion constant for the species u
  • Diffusion constant for the species v
  • Number of time-steps
  • Physical constant K
  • Physical constant F
openfpm_init(&argc,&argv);
// domain
Box<2,double> domain({0.0,0.0},{2.5,2.5});
// grid size
size_t sz[2] = {128,128};
// Define periodicity of the grid
periodicity<2> bc = {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 = 15000;
// K and F (Physical constant in the equation)
double K = 0.055;
double F = 0.03;
[v_transform metafunction]
Boundary conditions.
Definition common.hpp:22

Here we create 2 distributed grid in 2D Old and New. In particular because we want that the second grid is distributed across processors in the same way we pass the decomposition of the Old grid to the New one in the constructor with Old.getDecomposition(). Doing this, we force the two grid to have the same decomposition.

// New grid with the decomposition of the old grid
grid_dist_id<2, double, aggregate<double,double>> New(Old.getDecomposition(),sz,g);
// spacing of the grid on x and y
double spacing[2] = {Old.spacing(0),Old.spacing(1)};

We use the function init to initialize U and V on the grid Old

init(Old,New,domain);

Time stepping

After initialization, we first synchronize the ghost part of the species U and V for the grid, that we are going to read (Old). In the next we are going to do 15000 time steps using Eulerian integration

Because the update step of the Laplacian operator from \( \frac{\partial u}{\partial t} = \nabla u + ... \) discretized with eulerian time-stepping look like

\( \delta U_{next}(x,y) = \delta t D_u (\frac{U(x+1,y) - 2U(x,y) + U(x-1,y)}{(\delta x)^2} + \frac{U(x,y+1) - 2U(x,y) + U(x,y-1)}{(\delta y)^2}) + ... \)

If \( \delta x = \delta y \) we can simplify with

\( U_{next}(x,y) = \frac{\delta t D_u}{(\delta x)^2} (U(x+1,y) + U(x-1,y) + U(x,y-1) + U(x,y+1) -4U(x,y)) + ... \) (Eq 2)

The specie V follow the same concept while for the \( ... \) it simply expand into

\( - \delta t uv^2 - \delta t F(U - 1.0) \)

and V the same concept

See also
Ghost
Loop over grid points
Grid coordinates
VTK and visualization
// 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]);
for (size_t i = 0; i < timeSteps; ++i)
{
auto it = Old.getDomainIterator();
while (it.isNext())
{
auto key = it.get();
// update based on Eq 2
New.get<U>(key) = Old.get<U>(key) + uFactor * (
Old.get<U>(key.move(x,1)) +
Old.get<U>(key.move(x,-1)) +
Old.get<U>(key.move(y,1)) +
Old.get<U>(key.move(y,-1)) +
-4.0*Old.get<U>(key)) +
- deltaT * Old.get<U>(key) * Old.get<V>(key) * Old.get<V>(key) +
- deltaT * F * (Old.get<U>(key) - 1.0);
// update based on Eq 2
New.get<V>(key) = Old.get<V>(key) + vFactor * (
Old.get<V>(key.move(x,1)) +
Old.get<V>(key.move(x,-1)) +
Old.get<V>(key.move(y,1)) +
Old.get<V>(key.move(y,-1)) -
4*Old.get<V>(key)) +
deltaT * Old.get<U>(key) * Old.get<V>(key) * Old.get<V>(key) +
- deltaT * (F+K) * Old.get<V>(key);
// 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 100 time step we output the configuration for
// visualization
if (i % 100 == 0)
{
Old.ghost_get<U,V>();
Old.write_frame("output",count);
count++;
}
}

Finalize

Deinitialize the library

openfpm_finalize();