-- FreeFem++ v4.6 (Thu Apr  2 15:47:38 CEST 2020 - git v4.6)
 Load: lg_fem lg_mesh lg_mesh3 eigenvalue 
    1 : //  chemotaxis.edp
    2 : //
    3 : //  Discussion:
    4 : //
    5 : //    Chemotaxis model for biological pattern generation.
    6 : //
    7 : //  Location:
    8 : //
    9 : //    http://people.sc.fsu.edu/~jburkardt/freefem++/chemotaxis/chemotaxis.edp
   10 : //
   11 : //  Modified:
   12 : //
   13 : //    16 May 2015
   14 : //
   15 : //  Author:
   16 : //
   17 : //    Florian De Vuyst
   18 : //
   19 : //  Reference:
   20 : //
   21 : //    Florian De Vuyst,
   22 : //    Numerical modeling of transport problems using freefem++ software -
   23 : //    with examples in biology, CFD, traffic flow and energy transfer,
   24 : //    HAL id: cel-00842234
   25 : //    https://cel.archives-ouvertes.fr/cel-00842234
   26 : //
   27 : 
   28 : //
   29 : // Parameter definition
   30 : //
   31 : real r = 38.05;
   32 : real alpha = 285;
   33 : real rhoinf =1;
   34 : real D = 0.25;
   35 : real s = 1;
   36 : 
   37 : real dt = 0.05;
   38 : //
   39 : //  Define the mesh and plot it.
   40 : //
   41 : real Lx = 3.5;
   42 : real Ly = 4;
   43 : 
   44 : mesh Th=square(40, 40, [Lx*x, Ly*y]);
   45 : 
   46 : plot ( Th, wait = 0, ps = "chemotaxis_mesh.ps" );
   47 : //
   48 : //  Define the finite element spaces.
   49 : //
   50 : fespace Vh(Th, P1);
   51 : Vh rho;
   52 : Vh rhoold;
   53 : Vh rhotest;
   54 : Vh u1;
   55 : Vh u2;
   56 : Vh output;
   57 : Vh sigma;
   58 : Vh sigmaold;
   59 : Vh sigmatest;
   60 : 
   61 : fespace Wh(Th, P2);
   62 : Wh c;
   63 : Wh cold;
   64 : Wh ctest;
   65 : //
   66 : //  Start from the unstable constant state.
   67 : //
   68 : rho = rhoinf;
   69 : rhoold = rho;
   70 : c = rhoinf/(1+rhoinf);
   71 : cold=c;
   72 : u1 = alpha * dx(c);
   73 : u2 = alpha * dy(c);
   74 : //
   75 : //  Define the weak form of the equation for the density.
   76 : //
   77 : problem step1(rho, rhotest) =
   78 :   int2d(Th)(rho*rhotest /dt)
   79 : - int2d(Th)( convect([u1,u2], -dt, rhoold)*rhotest /dt)
   80 : + int2d(Th)(D*dx(rho)*dx(rhotest)+D*dy(rho)*dy(rhotest))
   81 : - int2d(Th)(r*s*rhoold*rhoinf*rhotest)
   82 : + int2d(Th)(r*s*rhoold*rho*rhotest);
   83 : //
   84 : //  Define the weak form of the equation for sigma and c.
   85 : //
   86 : problem step2([sigma, c], [sigmatest, ctest]) =
   87 :   int2d(Th)( sigma*sigmatest /dt)
   88 : - int2d(Th)( sigmaold*sigmatest /dt)
   89 : - int2d(Th)(alpha*dx(c)*dx(sigmatest)+alpha*dy(c)*dy(sigmatest))
   90 : + int2d(Th)(c*ctest /dt)
   91 : - int2d(Th)(cold*ctest /dt)
   92 : + int2d(Th)(dx(c)*dx(ctest) + dy(c)*dy(ctest))
   93 : - int2d(Th)(s*rho/(1+rho)*ctest)
   94 : + int2d(Th)(s*c*ctest);
   95 : 
   96 : for ( int it = 0; it < 200; it++ ) 
   97 : {
   98 :   cout << "it = " << it << endl;
   99 : //
  100 : // Step 1
  101 : //
  102 :   step1;
  103 :   rhoold = rho;
  104 :   sigmaold = log(rhoold);
  105 : //
  106 : //  Step 2
  107 : //
  108 :   step2;
  109 :   u1 = alpha * dx(c);
  110 :   u2 = alpha * dy(c);
  111 :   rho = exp(sigma);
  112 :   rhoold = rho;
  113 :   cold = c;
  114 : 
  115 :   if ( ( it % 10 ) == 0 && 0 )
  116 :   {
  117 :     plot ( c, nbiso = 60, value=1, wait=1 );
  118 :   }
  119 : 
  120 :   if ( ( it % 10 ) == 0 )
  121 :   {
  122 :     plot ( rho, nbiso = 40, grey = 0, fill = 1, value = 1, wait = 0 );
  123 :   }
  124 : }
  125 : 
  126 : cout << "cmin = " << c[].min << " cmax = " << c[].max << endl;
  127 : cout << "rhomin = " << rho[].min << " rhomax = " << rho[].max << endl;
  128 : //
  129 : plot( rho, nbiso=60, grey=1, fill=1, value=0, wait=0, ps="chemotaxis_rho.ps");
  130 : plot( c , nbiso=60, grey=1, fill=1, value=0, wait=0, ps="chemotaxis_c.ps");
  131 : //
  132 : //  Terminate.
  133 : //
  134 : cout << "\n";
  135 : cout << "CHEMOTAXIS:\n";
  136 : cout << "  Normal end of execution.\n";
  137 : 
  138 :  sizestack + 1024 =8560  ( 7536 )

  -- Square mesh : nb vertices  =1681 ,  nb triangles = 3200 ,  nb boundary edges 160
it = 0
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.09637e-11  max 8.47711e-12
          min 0.5  max 0.5
it = 1
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -7.85859e-12  max 9.1373e-12
          min 0.5  max 0.5
it = 2
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -9.6282e-12  max 1.04171e-11
          min 0.5  max 0.5
it = 3
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -9.70923e-12  max 1.1366e-11
          min 0.5  max 0.5
it = 4
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.00795e-11  max 1.42071e-11
          min 0.5  max 0.5
it = 5
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.05304e-11  max 1.55919e-11
          min 0.5  max 0.5
it = 6
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.14817e-11  max 1.70487e-11
          min 0.5  max 0.5
it = 7
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.31803e-11  max 2.2997e-11
          min 0.5  max 0.5
it = 8
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.63985e-11  max 2.58441e-11
          min 0.5  max 0.5
it = 9
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -2.02408e-11  max 3.20358e-11
          min 0.5  max 0.5
it = 10
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -2.27881e-11  max 3.79747e-11
          min 0.5  max 0.5
it = 11
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -2.68275e-11  max 4.68794e-11
          min 0.5  max 0.5
it = 12
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -3.34701e-11  max 5.68125e-11
          min 0.5  max 0.5
it = 13
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -3.95212e-11  max 6.85479e-11
          min 0.5  max 0.5
it = 14
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -4.53109e-11  max 8.55515e-11
          min 0.5  max 0.5
it = 15
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -5.73768e-11  max 1.03369e-10
          min 0.5  max 0.5
it = 16
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -6.79968e-11  max 1.2416e-10
          min 0.5  max 0.5
it = 17
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -8.24243e-11  max 1.50304e-10
          min 0.5  max 0.5
it = 18
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.01557e-10  max 1.80087e-10
          min 0.5  max 0.5
it = 19
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.22358e-10  max 2.2174e-10
          min 0.5  max 0.5
it = 20
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.50496e-10  max 2.64234e-10
          min 0.5  max 0.5
it = 21
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.85798e-10  max 3.23317e-10
          min 0.5  max 0.5
it = 22
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -2.27496e-10  max 3.99089e-10
          min 0.5  max 0.5
it = 23
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -2.81776e-10  max 4.83382e-10
          min 0.5  max 0.5
it = 24
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -3.45312e-10  max 5.91505e-10
          min 0.5  max 0.5
it = 25
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -4.26765e-10  max 7.21045e-10
          min 0.5  max 0.5
it = 26
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -5.22856e-10  max 8.78994e-10
          min 0.5  max 0.5
it = 27
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -6.43224e-10  max 1.07217e-09
          min 0.5  max 0.5
it = 28
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -7.91816e-10  max 1.30896e-09
          min 0.5  max 0.5
it = 29
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -9.71631e-10  max 1.5989e-09
          min 0.5  max 0.5
it = 30
  -- Solve : 
          min 1  max 1
  -- Solve : 
          min -1.19708e-09  max 1.94759e-09
          min 0.5  max 0.5