Upload of simple mPDeC code

DeC_Patankar_Notebook/DeC_procedure.jl

+using LinearAlgebra
+
+
+include("Lagrange_polynomials.jl")
+
+
+function dec_correction( IC,  M_sub::Int, K_corr::Int, y_0, ts)
+
+    #including production destruction functions
+    if IC=="linear"
+      include("Linear_model.jl")
+    elseif IC=="robertson"
+      include("Robertson_model.jl")
+    elseif IC=="non_linear"
+      include("Nonlinear_model.jl")
+    else
+      println("Not valid IC")
+    end
+
+    #Initial data, time(0), timesteps, dimension of the system
+    t_0 = ts[1] 
+    N_time = length(ts)-1
+    dim=length(y_0)
+    #Variables, U is the solution for all timesteps
+    U=zeros(dim, N_time+1)
+    #Variables of the DeC procedure: one for each subtimesteps and one for previous correction up, one for current correction ua
+    u_p=zeros(dim, M_sub+1)
+    u_a=zeros(dim, M_sub+1)
+
+    #Matrices of production and destructions applied to up at different subtimesteps
+    prod_p=zeros(dim,dim,M_sub+1)
+    dest_p=zeros(dim,dim,M_sub+1)
+
+    #coefficients for time integration
+    Theta=compute_theta(M_sub)
+
+    U[:,1]=y_0
+    delta_t=zeros(N_time)
+
+    #Timesteps loop
+    for it=2: N_time+1
+        #computing actual delta t
+        delta_t[it-1]=ts[it]-ts[it-1]
+
+        # Subtimesteps loop to update variables at iteration 0
+        for m=1:M_sub+1
+         u_a[:,m]=U[:,it-1]
+         u_p[:,m]=U[:,it-1]
+        end
+        #Loop for iterations K
+        for k=2:K_corr+1
+            # Updating previous iteration variabls
+            u_p=copy(u_a)
+            #Loop on subtimesteps to compute the production destruction terms
+            for r=1:M_sub+1
+                prod_p[:,:,r], dest_p[:,:,r]=prod_dest(u_p[:,r])
+            end
+            # Loop on subtimesteps to compute the new variables
+            for m=2:M_sub+1
+                u_a[:,m]=patanker_type_dec(prod_p, dest_p, delta_t[it-1], m, M_sub, Theta, u_p, dim)
+            end
+        end
+        # Updating variable solution
+        U[:,it]=u_a[:,M_sub+1]
+    end
+    return delta_t, U
+end
+
+
+"""prod_p and dest_p dim x dim
+delta_t is timestep length
+m_substep is the subtimesteps for which we are solving the system
+M_sub is the maximum number of subtimesteps (0,..., M)
+theta coefficient vector matrix  in  M_sub+1 x M_sub+1, first index is the index we are looping on,
+the second is the one we are solving for, theta_r,m= int_t0^tm phi_r
+u_p is the previous correction solution for all subtimesteps in dim x (M_sub+1) """
+function patanker_type_dec(prod_p, dest_p, delta_t,  m_substep::Int, M_sub, theta, u_p, dim)
+#This function builds and solves the system u^{(k+1)}=Mu^{(k)} for a subtimestep m
+
+mass=Matrix{Float64}(I, dim, dim);
+#println(mass)
+for i=1:dim
+    for r=1: M_sub+1
+        if theta[r,m_substep]>0
+            for j= 1: dim
+              mass[i,j]=mass[i,j]- delta_t*theta[r,m_substep]*(prod_p[i,j,r]/u_p[j,m_substep])
+              mass[i,i]=mass[i,i]+ delta_t*theta[r,m_substep]*(dest_p[i,j,r]/u_p[i,m_substep])
+            end
+        elseif theta[r,m_substep]<0
+            for j= 1: dim
+              mass[i,i]=mass[i,i]- delta_t*theta[r,m_substep]*(prod_p[i,j,r]/u_p[i,m_substep])
+              mass[i,j]=mass[i,j]+ delta_t*theta[r,m_substep]*(dest_p[i,j,r]/u_p[j,m_substep])
+            end
+        end
+    end
+end
+return mass\u_p[:,1]
+end
+

DeC_Patankar_Notebook/Lagrange_polynomials.jl

+using SymPy
+
+
+function generate_poly(s, nodes,deg)
+#Generates Lagrangian polynomial in symbolic
+
+aux=[]
+
+for k=1:deg+1
+    append!(aux, 1);
+    for j=1:deg+1
+        if j != k
+            aux[k]=aux[k].*(s-nodes[j]);
+        end
+    end
+end
+for j=1:deg+1
+    aux[j]= aux[j]/SymPy.N(subs(aux[j],s,nodes[j]));
+end
+return aux
+end
+
+function compute_theta(deg)
+# Compute the integral of lagrange polynomials from 0 to interpolation point theta[r,m]= \int_{t^0}^{t^m} \phi_r
+s=Sym("s");
+nodes=range(0,stop=1,length=deg+1);
+
+poly=generate_poly(s,nodes,deg)
+
+theta=zeros(deg+1, deg+1)
+
+for r=1:deg+1
+    for m=1:deg+1
+        theta[r,m]=integrate(poly[r],(s,0,nodes[m]))
+    end
+end
+
+return theta
+
+end

DeC_Patankar_Notebook/Linear_model.jl

+function prod_dest(c)
+#Production destruction matrices for a linear model
+a=parameters()
+d= length(c)
+p=zeros(d,d)
+d=zeros(d,d)
+p[1,2]=c[2]
+d[2,1]=c[2]
+p[2,1]=a*c[1]
+d[1,2]=a*c[1]
+return p, d
+end
+
+"""Exact Solution for the linear Problem """
+function exact_solution(c0, time)
+a= parameters()
+c1_inf= sum(c0)/(a+1)
+c = c0[1]/c1_inf -1.0
+
+c1 = (1+c*exp(-(a+1)*time))*c1_inf
+c2=sum(c0)-c1
+return [c1,c2]
+
+end
+
+
+function parameters()
+a=5.0
+return a
+end

DeC_Patankar_Notebook/Nonlinear_model.jl

+function prod_dest(c)
+#Production destruction matrices for a nonlinear model
+a=parameters_non()
+d= length(c)
+p=zeros(d,d)
+d=zeros(d,d)
+d[1,2]=c[1]*c[2]/(c[1]+1.)
+p[2,1]=c[1]*c[2]/(c[1]+1.)
+p[3,2]=a*c[2]
+d[2,3]=a*c[2]
+return p, d
+end
+
+function parameters_non()
+a=0.3
+return a
+end

DeC_Patankar_Notebook/Robertson_model.jl

+function prod_dest(c)
+#Production destruction matrices for a Robertson model
+d= length(c)
+p=zeros(d,d)
+d=zeros(d,d)
+p[1,2]=10^4*c[2]*c[3]
+d[2,1]=10^4*c[2]*c[3]
+p[2,1]=0.04*c[1]
+d[1,2]=0.04*c[1]
+p[3,2]=3*10^7*c[2]^2
+d[2,3]=3*10^7*c[2]^2
+return p, d
+end
+

DeC_Patankar_Notebook/Time_integration.jl

+using Plots
+using LinearAlgebra
+
+include("Lagrange_polynomials.jl")
+include("DeC_procedure.jl")
+
+#3 possible tests
+test="robertson" #     "non_linear"#  "linear"  #
+
+#number of timesteps
+N_time = 100
+# order d => subtimesteps M=order-1, corrections = order
+order=4
+M_sub=order-1
+K_iter=order
+
+# Test data: initial conditions, times
+if test=="linear"
+  include("Linear_model.jl")
+  u0 = [0.9;0.1]
+  T_0=0
+  T_final = 1.75
+  ts=range(T_0,stop=T_final,length=N_time+1)
+elseif test=="robertson"
+  epsilon= Base.eps()
+  u0=[1.0-2.0*epsilon;epsilon;epsilon]
+
+  T_0 = 10^-6
+  T_final = 10^10
+
+  "Variable timesteps in logaritmic scale"
+  logdt0 = (log(T_final)-log(T_0))/N_time
+  dt0=exp(logdt0)
+  ts= zeros(N_time+1)
+  ts[1]=T_0
+  for it=1:N_time
+    ts[it+1]=T_0*dt0^it
+  end
+
+  include("Robertson_model.jl")
+elseif test=="non_linear"
+  u0 = [9.98;0.01;0.01]
+  T_0=0
+  T_final = 30.
+  ts=range(T_0,stop=T_final,length=N_time+1)
+  include("Nonlinear_model.jl")
+else
+  println("Not valid IC")
+end
+
+#Computing the dec patankar method 
+dts, U = dec_correction( test,  M_sub, K_iter, u0, ts ) 
+
+#Plotting the solutions
+if test=="robertson"
+#For robertson test case, we multiply the second variable times 10^4 to be visible
+  U[2,:]=10^4*U[2,:]
+  plot(ts,transpose(U), xaxis=:log)
+else
+  plot(ts,transpose(U))
+end  
+plot!(title="Test $(test)",
+     xaxis="Time",yaxis="Values", legend=:best)
+savefig("../Figures/$(test)_N$(N_time).pdf")