//  lotka_volterra.edp
//
//  Discussion:
//
//    Solve the reaction-diffusion equation of the Lotka-Volterra 
//    predator-prey model.
//
//    The differential reaction system is solved by a Strang splitting.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/lotka_volterra/lotka_volterra.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 May 2020
//
//  Author:
//
//    Florian De Vuyst
//
//  Reference:
//
//    Florian De Vuyst,
//    Numerical modeling of transport problems using freefem++ software -
//    with examples in biology, CFD, traffic flow and energy transfer,
//    HAL id: cel-00842234
//    https://cel.archives-ouvertes.fr/cel-00842234
//
cout << "\n";
cout << "lotka_volterra\n";
cout << "  FreeFem++ version\n";
cout << "  Simulate the time evolution of the spatial distribution of\n";
cout << "  predator and prey species in a 2D domain.\n";
//
//  Set parameters.
//
real alpha = 0.5;
real t = 0.0;
real dt = 0.2;
real nu = 0.001;
real mu = 0.0001;
real theta = 0.49;
int it = 0;
//
//  Th: the mesh.
//
mesh Th = square ( 40, 40 );
//
//  Plot the mesh.
//
plot ( Th, ps = "lotka_volterra_mesh.ps" );
//
//  Define the finite element spaces.
//
fespace Vh ( Th, P2 );
fespace Wh ( Th, P1 );
Vh u, uh, uold;
Vh v, vh, vold;
Wh dxu, dyu, dxv, dyv;
//
//  Define the weak form of the U equation.
//
problem heatu (u, uh, init=it) =
  int2d (Th) (u*uh/dt)
- int2d (Th) (uold*uh/dt)
+ int2d (Th) (nu*(1.0-theta)*dx(u)*dx(uh) + nu*(1.0-theta)*dy(u)*dy(uh))
+ int2d (Th) (nu*(theta)*dx(uold)*dx(uh) + nu*(theta)*dy(uold)*dy(uh));
//
//  Define the weak form of the V equation.
//
problem heatv (v, vh, init=it) =
  int2d (Th) (v*vh/dt)
- int2d (Th) (vold*vh/dt)
+ int2d (Th) (mu*(1.0-theta)*dx(v)*dx(vh)+mu*(1.0-theta)*dy(v)*dy(vh))
+ int2d (Th) (mu*(theta)*dx(vold)*dx(vh)+mu*(theta)*dy(vold)*dy(vh));
//
// Initializing
//
u = 1.0 + 0.5 * cos ( 3.0 * pi * x ) * sin ( 3.0 * pi * y );
uold = u;
v = 1.0 + 0.5 * sin ( 3.0 * pi * x ) * cos ( 3.0 * pi * x );
vold = v;
//
//  Big loop in time
//  Fractional step method
//
for ( it = 0; it < 200; it++ ) 
{
  for ( int subit = 0; subit < 1; subit++ )
  {
    t = t + dt;
//
//  Step a. Solve the reaction system on ( dt /2)
//
    u = u * exp ( alpha*(v-1.0)*dt*0.25 );
    v = v * exp ( (1.0-u)*dt*0.5 );
    u = u * exp ( alpha*(v-1.0)*dt*0.25 );

    uold = u; 
    vold = v;
//
// Step b. Solve the diffusion system on ( dt )
//
    heatu; 
    uold = u;
    heatv; 
    vold = v;
//
// Step c. Solve the reaction system on ( dt /2)
//
    u = u * exp ( alpha * ( v - 1.0 ) * dt * 0.25 );
    v = v * exp ( ( 1.0 - u ) * dt * 0.5 );
    u = u * exp ( alpha * ( v - 1.0 ) * dt * 0.25 );
    uold = u; 
    vold = v;
  }
//
//  Display current solution to the interactive screen.
//
  plot ( u, nbiso = 40, value = true, fill = true, cmm = "U species" );
  plot ( v, nbiso = 40, value = true, fill = true, cmm = "V species" );
}
//
//  Save the final images.
//
  plot ( u, nbiso = 40, value = true, fill = true, 
    ps = "lotka_volterra_u.ps", cmm = "U species" );

  plot ( v, nbiso = 40, value = true, fill = true, 
    ps = "lotka_volterra_v.ps", cmm = "V species" );
//
//  Terminate.
//
cout << "\n";
cout << "lotka_volterra:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );