IMPETUS – Interactive MultiPhysics Simulation Environment
Tutorials
Tutorial 0: Using the Tutorial
Tutorial 1: Creating a simple Molecular Dynamics Simulation
Tutorial 2: Solving the Diffusion equation using Interactive Field
Tutorial 3: Simple Chemotaxis Model
Full PDF tutorial is also available.
Download IMPETUS and Tutorials.
Tutorial 3: Simple Chemotaxis Model
3.1 Introduction
This section of the tutorial will guide the users to couple the “CellSpace” MID with the “InteractiveField” MID into one simulation. We will demonstrate a numerical model that simulates chemotaxis. Particles from “CellSpace” will detect the concentration field generated by the “InteractiveField” and migrate along the maximum concentration gradient.
This section of the tutorial is located in:
Tutorials/tutorialChemotaxis/
Include the header for the simple molecular dynamics simulation by inserting the following line in main.cpp
#include "../Tutorials/tutorialChemotaxis/chemotaxis.h"
3.2 More InteractiveField Functions
·
To get the concentration 
at a certain coordinate x0, y0, z0,the user should use:
double concentration = c -> getConcentration(x0, y0, z0);
·
To get the gradient of the concentration 
at a certain coordinate x0, y0, z0, the user should use:
vec3 gradient = c -> getGradient(x0, y0, z0);
3.3 Building a Cross-interactive Function using “getPartList”:
The function “getPartList” will use the coordinate of the particles in the “CellSpace” to acquire the gradient vector of the diffusion. The particles will follow this gradient to migrate towards the highest concentration point.
![Text Box: class Migration_tutorial : public getPartList {
public:
double threshold;
double f;
Migration_tutorial(vamde::CellSpace *_s, vamde::InteractiveField * _c) : getPartList(_s) , c(_c){
threshold = 1 ;
f = 0.25;
}
private:
vamde::InteractiveField * c;
void action(Particle *i){
double concentration = c -> getConcentration(i->x[0],i->x[1],i->x[2]);
vec3 gradient = c -> getGradient(i->x[0],i->x[1],i->x[2]);
double dr[3];
double rr = 0;
for (int d=0; d<DIM; d++) {
dr[d] = gradient.r[d];
rr += dr[d] * dr[d];
}
double x =sqrt(rr);
double drhat[3];
for (int d=0; d<DIM; d++) {
drhat[d] = dr[d] / x;
}
if (x > threshold) {
for (int d=0; d<DIM; d++){
i->F[d] += f * drhat[d];
}
}
}
};](tutword3_files/image003.gif)
· The function is called as:
Migration_tutorial mg0(s0, c0);
· It is used as:
mg0.apply();
3.4 Creating Sinks and Sources
In this example, we introduce a sink and a source in the concentration field. The source is simulated by constantly setting the concentration at the point where the source is located at a specific value. Here, we use the value of 10. The sink is simulated by constantly setting the concentration to -10). The two functions are placed in a “while” loop.
double x_mid = (param->gridinfo->world.lo[0] + param->gridinfo->world.hi[0]) /2;
double y_mid = (param->gridinfo->world.lo[1] + param->gridinfo->world.hi[1]) /2;
double z_mid = (param->gridinfo->world.lo[2] + param->gridinfo->world.hi[2]) /2;
c0->setConcentration(x_mid+0.4*x_mid,y_mid,z_mid,10 );
c0->setConcentration(x_mid-0.4*x_mid,y_mid,z_mid,-10 );
By combining everything that we discussed so far, we can create the following chemotaxis model:
![Text Box: void runSimulation(vamde::Parameters * param) {
vamde::TimeKeeper * tk;
tk = new vamde::TimeKeeper(param);
/// parameters:
double delta_t = param-> delta_t;
double end_step = param-> end_step;
double & t = param->clock->t;
int & step = param->clock->step;
/// Initiate clock
param->clock->set_step(0);
vamde::CellSpace ** s;
s = new vamde::CellSpace * [param->n_cell];
s[0] = new vamde::CellSpace("Tutorials/tutorial3chemotaxis/input/cell0.input",param);
s[0]->cfgwriter->set_atom_name( "Cs");
/// Reset the unique ID of all particles
for (int n=0; n<param->n_cell; n++) {
s[n]->refreshCell();
}
vamde::InteractiveField ** c;
c = new vamde::InteractiveField * [param->n_cont];
vamde::InteractiveField::Cinit continit;
continit.readinput("Tutorials/tutorial3chemotaxis/input/cont0.input");
c[0] = new vamde::InteractiveField(continit,param);
c[0]->diffusivity = 0.1;
double x_mid = (param->gridinfo->world.lo[0] + param->gridinfo->world.hi[0]) /2;
double y_mid = (param->gridinfo->world.lo[1] + param->gridinfo->world.hi[1]) /2;
double z_mid = (param->gridinfo->world.lo[2] + param->gridinfo->world.hi[2]) /2;
c[0]->setConcentration(x_mid+0.4*x_mid,y_mid,z_mid, 10 );
c[0]->setConcentration(x_mid-0.4*x_mid,y_mid,z_mid,-10 );
c[0]->setAllGlobalBoundarySink();
for (int n=0; n<param->n_cont; n++) {
c[n]->cfgwriter->print();
}
LenardJones3 lenardjones(s,param->n_cell );
ZeroForce3 zeroforce(s,param->n_cell);
LeapFrog3 leapfrog(s,param->n_cell);
Viscosity3 viscosity(s,param->n_cell);
for (int n=0; n<param->n_cell; n++) {
s[n]->deleteBorders();
s[n]->cfgwriter->print();
}
while (step < end_step) {
/// Print Runtime progress
tk->print_progress(10);
/// Advace the clock by one step
param->clock->advance();
/// LeapFrog step 1
leapfrog.step(1);
for (int n=0; n<param->n_cell; n++) {
/// delete ghost particles
s[n]->deleteBorders();
/// move particles to new cells after integration before moving interprocessor
s[n]->move_particles_to_the_correct_cells();
/// move particles across processors.
s[n]->moveParticles();
/// copy new ghost particles
s[n]->copyParticles();
}
/// Set All forces to zero
zeroforce.loopLocalParticles();
/// Apply Lenard Jones Pair potential
lenardjones.loopNeighbors();
viscosity.loopLocalParticles();
/// LeapFrog step 2
leapfrog.step(2);
/// print Cfgs
for (int n=0; n<param->n_cell; n++) {
s[n]->cfgwriter->print();
}
for (int n=0; n<param->n_cont; n++) {
c[n]->copyParticles();
}
for (int n=0; n<param->n_cont; n++) {
c[n]->Iterate();
}
/// continuum
c[0]->cfgwriter->print();
c[0]->setConcentration(x_mid+0.4*x_mid,y_mid,z_mid, 10 );
c[0]->setConcentration(x_mid-0.4*x_mid,y_mid,z_mid, -10 );
Migration_tutorial mg0(s[0], c[0]);
mg0.apply();
}
}](tutword3_files/image004.gif)
Results are shown in EXAMPLES.