How to speed up NDSolve when there are ~400k equations?

Question

I'm trying to use NDSolve to solve a set of ~400k equations. The problem I met is: when I ran the evaluation, Mathematica (version 11.0) raised an error:

NDSolve::ntdv: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"}.

However, after I added the option Method->{"EquationSimplification"->"Residual"} into the evaluation, it took almost 24h to get the anwsers.

So, what does the option Method->{"EquationSimplification"->"Residual"} really means? Are there other ways I can get the equations solved correctly and quickly, maybe in 1h?

Below are my simlified codes:

(*Equations*)
EqD[nC_,nN_] := c[nC,nN]'[t];
Eqk[nC_,nN_] := Sum[b[nC-i,nN-(1-i),i,1-i]*c[nC-i,nN-(1-i)][t]-(b[nC,nN,i,1-i]+a[nC,nN,i,1-i])*c[nC,nN][t]+a[nC+i,nN+(1-i),i,1-i]*c[nC+i,nN+(1-i)][t],{i,{0,1}}];
Eq[nC_,nN_] := EqD[nC,nN] == Eqk[nC,nN];
(*Coefficients*)
b[nC_,nN_,i_,j_] := 4\[Pi]*(R[nC,nN]+R[i,j])/\[CapitalOmega]*D\[Alpha][i,j]*c[i,j][t];
a[nC_,nN_,i_,j_] := 4\[Pi]*(R[nC-i,nN-j]+R[i,j])/\[CapitalOmega]*D\[Alpha][i,j]*Exp[-((F[nC-i,nN-j]+F[i,j]-F[nC,nN])/(k*T))]
(*Other definitions needed*)
R[nC_,nN_] := (3(nC+nN)*\[CapitalOmega]/(4\[Pi]))^(1/3)
D\[Alpha][i_,j_] := 0;D\[Alpha][1,0] := 8.0*10^-4;D\[Alpha][0,1] := 1.0*10^-4;
F[nC_,nN_] := -k*T*Log[(nC+nN)!/(nC!*nN!)]+(nC*Fc[nC+nN]+nN*Fn[nC+nN])/(nC+nN);
(*\[Pi],\[CapitalOmega],k and T are constants;Fc[nC_] and Fn[nN_] are two functions related to R[nC,0] and R[0,nN]*)

Are there any useful suggestions to speed up my evaluation? Thanks a lot!!!

Answer 1

I tested the code as it is without optimization. I give a report and general recommendations: Test 1

In[11]:= Equations[tMin, tMax, 200, 20]; // AbsoluteTiming

Out[11]= {7.76804, Null}

Test 2

In[12]:= Equations[tMin, tMax, 200, 200]; // AbsoluteTiming

Out[12]= {86.8025, Null}

Test 3

In[13]:= Equations[tMin, tMax, 2000, 20]; // AbsoluteTiming

Out[13]= {109.039, Null}

Test 4

In[14]:= Equations[tMin, tMax, 2000, 200]; // AbsoluteTiming

Out[14]= {2174.46, Null}

And so it took 36 minutes, 30-50% of the CPU and 12GB of memory (maximum). Most of the time 1 kernel worked. I think that we can reduce the time to 3-5 minutes. Version report:

 $Version

Out[]= "12.0.0 for Microsoft Windows (64-bit) (April 6, 2019)"

Let's start the optimization. Here we are dealing with the Riccati matrix equation $$c'=c.A.c+B1.c+c.B2...(1) $$ where c is a 2000x200 matrix. Therefore we used Method -> {"EquationSimplification" -> "Residual"to solve it. Then it takes about 36 min on my ASUS laptop. We can reduce time by simple replacement t->tmax t. The corresponding code is

PrepareEquations[tMin_, tMax_, NoC_, NoN_] := 
  Module[{i, 
    j,(*C and N species*)\[Alpha]C = 1,(*C monomer*)\[Alpha]N = 
     1,(*N monomer*)nCmin = 0, nCmax = NoC,(*Number of C atoms*)
    nNmin = 0, nNmax = NoN,(*Number of N atoms*)
    c0 = 10^-30(*For numerical stability,avoid null value*)}, 
   Print["Function: PrepareEquations"];
   Clear[c, a, b, F, Fc, Fn];
   Cinit = 0.00439673;(*Initial C concentration*)
   Ninit = 0.00761477;(*Initial N concentration*)
   k = 1.380649*10^-23;(*Boltzmann constant,unit J/K*)
   T = 723;(*Temperature,unit K*)
   R[nC_, nN_] := ((3*(nC + nN)*0.0118)/(4 \[Pi]))^(1/3);(*Radius*)
   D\[Alpha][i_, j_] := 0;(*Diffusion coefficients*)
   D\[Alpha][\[Alpha]C, 0] := 8.0*10^-4;
   D\[Alpha][0, \[Alpha]N] := 1.0*10^-4;
   k\[Alpha][nC_, nN_, i_, j_] := 
    4 \[Pi]*(R[nC, nN] + R[i, j])/0.0118*D\[Alpha][i, j];
   (*Calculation of energy*)RsC = 0.02780299;
   \[Sigma]C = 4.09324*10^-19;
   Fc[nC_] := 4 \[Pi]*(R[nC, 0])^2*\[Sigma]C*(1 - RsC/R[nC, 0]);
   \[Sigma]N = 5.54492*10^-20;
   RsN = -0.14758424;
   Fn[nN_] := 4 \[Pi]*(R[0, nN])^2*\[Sigma]N*(1 - RsN/R[0, nN]);
   F[nC_, 
     nN_] := -k*T*
      Log[(nC + nN)!/(nC!*nN!)] + (nC*Fc[nC + nN] + 
        nN*Fn[nC + nN])/(nC + nN);
   b[nC_, nN_, i_, j_] := 
    k\[Alpha][nC, nN, i, j]*c[i, j][t];(*Condensation rate*)
   a[nC_, nN_, i_, j_] := 
    k\[Alpha][nC - i, nN - j, i, j]*
     Exp[-((F[nC - i, nN - j] + F[i, j] - F[nC, nN])/(k*
           T))];(*Emission rate*)a[i_, j_, i_, j_] := 0;(*Monomer*)
   a[0, nN_, 1, j_] := 0;
   a[nC_, 0, i_, 1] := 0;
   Clear[CL, Eqk, EqD, Eq];
   With[{},(*Initial conditions*)CL[nC_, nN_] := c[nC, nN][tMin] == 0;
    (*Definition of the equilibrium equations*)
    EqD[nC_, nN_] := c[nC, nN]'[t];
    Eqk[nC_, nN_] := 
     Sum[b[nC - i, nN - (1 - i), i, 1 - i]*
        c[nC - i, nN - (1 - i)][
         t] - (b[nC, nN, i, 1 - i] + a[nC, nN, i, 1 - i])*
        c[nC, nN][t] + 
       a[nC + i, nN + (1 - i), i, 1 - i]*
        c[nC + i, nN + (1 - i)][t], {i, {0, 1}}];
    Eq[nC_, nN_] := EqD[nC, nN] == tMax Eqk[nC, nN];
    (*Definition of equilibrium concentrations for each special size*)
    Module[{nC, nN}, nC = nCmax;
     For[nN = 1, nN < nNmax, nN++, 
      Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]*
           c[nC, nN - 1][t] - (b[nC, nN, 0, 1] + a[nC, nN, 0, 1])*
           c[nC, nN][t] + 
          a[nC, nN + 1, 0, 1]*
           c[nC, nN + 1][t]) + (b[nC - 1, nN, 1, 0]*c[nC - 1, nN][t] -
           a[nC, nN, 1, 0]*c[nC, nN][t])];
     nN = nNmax;
     For[nC = 1, nC < nCmax, nC++, 
      Eqk[nC, nN] = (b[nC - 1, nN, 1, 0]*
           c[nC - 1, nN][t] - (b[nC, nN, 1, 0] + a[nC, nN, 1, 0])*
           c[nC, nN][t] + 
          a[nC + 1, nN, 1, 0]*
           c[nC + 1, nN][t]) + (b[nC, nN - 1, 0, 1]*c[nC, nN - 1][t] -
           a[nC, nN, 0, 1]*c[nC, nN][t])];
     nC = 0;
     For[nN = 2, nN < nNmax, nN++, 
      Eqk[nC, nN] = (b[nC, nN - 1, 0, 1]*
           c[nC, nN - 1][t] - (b[nC, nN, 0, 1] + a[nC, nN, 0, 1])*
           c[nC, nN][t] + 
          a[nC, nN + 1, 0, 1]*
           c[nC, nN + 1][t]) + (-b[nC, nN, 1, 0] c[nC, nN][t] + 
          a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t])];
     nN = 0;
     For[nC = 2, nC < nCmax, nC++, 
      Eqk[nC, nN] = (b[nC - 1, nN, 1, 0]*
           c[nC - 1, nN][t] - (b[nC, nN, 1, 0] + a[nC, nN, 1, 0])*
           c[nC, nN][t] + 
          a[nC + 1, nN, 1, 0]*
           c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1] c[nC, nN][t] + 
          a[nC, nN + 1, 0, 1]*c[nC, nN + 1][t])];
     nC = nCmax; nN = 0;
     Eqk[nC, 
       nN] = (-b[nC, nN, 0, 1]*c[nC, nN][t] + 
         a[nC, nN + 1, 0, 1]*c[nC, nN + 1][t]) + (b[nC - 1, nN, 1, 0]*
          c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t]);
     nC = 0; nN = nNmax;
     Eqk[nC, 
       nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + 
         a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (b[nC, nN - 1, 0, 1]*
          c[nC, nN - 1][t] - a[nC, nN, 0, 1]*c[nC, nN][t]);
     nC = nCmax; nN = nNmax;
     Eqk[nC, 
       nN] = (b[nC, nN - 1, 0, 1]*c[nC, nN - 1][t] - 
         a[nC, nN, 0, 1]*c[nC, nN][t]) + (b[nC - 1, nN, 1, 0]*
          c[nC - 1, nN][t] - a[nC, nN, 1, 0]*c[nC, nN][t]);
     nC = 0; nN = 0;
     CL[nC, nN] = c[nC, nN][tMin] == 0;
     Eqk[nC, nN] = 0;
     nC = \[Alpha]C; nN = 0;
     CL[nC, nN] = c[nC, nN][tMin] == Cinit + c0;
     Eqk[nC, 
       nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + 
         a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1]*
          c[nC, nN][t] + 
         a[nC, nN + 1, 0, 1]*
          c[nC, nN + 1][
           t]) - (Sum[(b[x, y, nC, nN] - a[x, y, nC, nN])*
           c[x, y][t], {x, 0, nCmax - 1}, {y, 0, nNmax}] + 
         Sum[(-a[nCmax, i, nC, nN]*c[nCmax, i][t]), {i, 0, nNmax}]);
     nC = 0; nN = \[Alpha]N;
     CL[nC, nN] = c[nC, nN][tMin] == Ninit + c0;
     Eqk[nC, 
       nN] = (-b[nC, nN, 1, 0]*c[nC, nN][t] + 
         a[nC + 1, nN, 1, 0]*c[nC + 1, nN][t]) + (-b[nC, nN, 0, 1]*
          c[nC, nN][t] + 
         a[nC, nN + 1, 0, 1]*
          c[nC, nN + 1][
           t]) - (Sum[(b[x, y, nC, nN] - a[x, y, nC, nN])*
           c[x, y][t], {x, 0, nCmax}, {y, 0, nNmax - 1}] + 
         Sum[(-a[i, nNmax, nC, nN]*c[i, nNmax][t]), {i, 0, nCmax}]);];
    Bc := 
     Flatten[ParallelTable[
        c[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]] // Evaluate];
   DistributeDefinitions[Eq, CL];
   Sys = Join[
       Flatten[ParallelTable[
         Eq[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]], 
       Flatten[ParallelTable[
         CL[nC, nN], {nC, 0, nCmax}, {nN, 0, nNmax}]]] // Flatten // 
     Evaluate;
   {Sys, Bc}];
tMin = 0; tMax = 7.75*10^5;

Here we have replaced only one equation Eq[nC_, nN_] := EqD[nC, nN] == tMax Eqk[nC, nN];.Now we call

PrepareEquations[tMin, tMax, 2000, 200]; // AbsoluteTiming

And it takes {536.24, Null}. Now we solve the system on {t, 0, 1}:

RSol = NDSolve[Sys, Bc, {t, tMin, 1}, 
    Method -> {"EquationSimplification" -> 
       "Residual"}]; // AbsoluteTiming

It takes only {1148.03, Null}. Thus, we saved 9 minutes with a simple replacement t->tmax t. However, if we can write the equation in matrix form, then we can use exponential Rosenbrock-type integrators. Then we can reduce the time to 1 minute.