kill(all)$

maxima:fvm:appendix:index notation and patankar(1980),convection-diffusion,upwind,version 1.0
peter.vlasschaert@gmail.com,31/01/2023

load packages :
-------------------

load(lrats);

definition : my product 
------------------------------
from :list to product

vt:[1, 2, 3, 4, 5];

 myprodlist(xs) := block([],
     prod: 1,
for i  in xs  do(
 prod : prod*i
),
return(prod/1)
)$

 myprodlist(vt);

use : factorial product  .

rem: product in python
------------------------------------
def list_product(numbers):
  product = 1
  for num in numbers:
    product *= num
  return product

numbers = [1, 2, 3, 4, 5]
result = list_product(numbers)
print(result)  # Output: 120

difference between : maxima you need block use return, maxma : ,python =,
                                   def name : python , name() :=block.

rem :

sa :[u[v1],Φ[v1],ρ[v1]] ,you can't use product form maxima.
→ product(sa[i],i,1,length(sa)) , give error

info: CV = Control Volume.
--------------------------------------
Patankar   : W,w,p,e;E
Index not. : i-1,i-1/2,i,i+1/2,i+1
---------------------------------------------------
distance : |wp| =Δx(i-1/2), |ep| = Δx(i+1/2)
distance : |Wp| =Δx(i-1),    |Ep| = Δx(i+1)
----------------------------------------------------
√ ∈ ,distance >0
-----------------------------------------------------

mathematical model : convection-diffusion

----------------------------------------------------------------

depends(Φ,x);

pp:m ='diff(Φ,x);

p1:'diff(u*ρ*Φ,x) = 'diff(D*m,x);

u = velocity [m/s]
ρ = density  [kg/m^3]
D = diffusion coef. [m^2/s]

-----------------------------------------------------------------

simp:false$

first step : integration over CV

p1a:integrate(lhs(p1),x,i-1/2,i+1/2) =integrate(rhs(p1),x,i-1/2,i+1/2); 

integration : CV , x=i-1/2  to x = i+1/2
-----------------------------------------------

→ first step :integration without boundary

p1b:(part(args(lhs(p1a))[1],1))=part(args(rhs(p1a))[1],1);

p1c1:makelist(part(lhs(p1b),j)[v1],j,length(lhs(p1b)));
p1c1a:ratsimp(myprodlist(p1c1));

p2a:at(p1c1a,v1=i+1/2) -at(p1c1a,v1=i-1/2);

p2a1::subst([v1=i+1/2],p1c1a) - subst([v1=i-1/2],p1c1a);

p1c2:makelist(part(rhs(p1b),j)[v1],j,length(rhs(p1b)));
p1c2a:ratsimp(myprodlist(p1c2));
p1c3a:ev(lratsubst([m[v1]=at('diff(Φ,x),x=v1)],p1c2a),nouns);

p2b:subst([v1=i+1/2],p1c3a) - subst([v1=i-1/2],p1c3a);

final result : integration "p1a"
-------------------------------------

p3:p2a=p2b;

some txt : Γ =D

use :F =ρ*u

let(u[i+1/2]*ρ[i+1/2],F[i+1/2]);
let(u[i-1/2]*ρ[i-1/2],F[i-1/2]);

p4:num(letsimp(lhs(p3)))=num(letsimp(rhs(p3)));

How to use : let  → letsimp
-------------------------------------
matchdeclare ([a], true)$
let (sin(a)^2, 1 - cos(a)^2);
letsimp (sin(x)^4);
outp :(cos (x))^4 - 2 cos (x) + 1
-------------------------------------
see :https://maxima.sourceforge.io/docs/manual/maxima_156.html#Item_003a-Rules_002fdeffn_002flet

disprule(all);

part(args(rhs(p4))[2],2);

values are : known,{ i-1,i,i+1}

diffusion part : finite difference 
---------------------------------------

find : differential  →

use : CDS = Central Difference Scheme

p5a:part(args(rhs(p4))[1],1,2) =(Φ[i]-Φ[i-1])/Δx(i-1);

p5b:part(args(rhs(p4))[2],2) =(Φ[i+1]-Φ[i])/Δx(i+1);

convection part : finite difference 
-----------------------------------------

p6a:Φ[i+1/2] =1/2*(Φ[i]+Φ[i+1] );
p6b:Φ[i-1/2] =1/2*(Φ[i]+Φ[i-1] );

combine :

p7:[p5a,p5b,p6a,p6b];

final : result 
---------------

p8:subst(p7,p4);

upwind scheme :  Φ[i+1/2], F =u*ρ
-------------------------------------------

u>0

if F[i+1/2] > 0 then Φ[i+1/2] :Φ[i] ;

u<0

if F[i+1/2] < 0 then Φ[i+1/2] :Φ[i+1] ;

patankar notation: (1980)

p9:F[i+1/2]*Φ[i+1/2] = Φ[i]*max(F[i+1/2],0) -Φ[i+1]*max(-F[i+1/2],0);

p9a:F[e]*Φ[e] =Φ[p]*max(F[i-1/2],0) -Φ[E]*max(-F[i-1/2],0);/*see figure 1*/

upwind scheme :  Φ[i-1/2],F=u*ρ
----------------------------------------

u>0

if F[i-1/2] > 0 then Φ[i-1/2] :Φ[i-1] ;

u<0

if F[i-1/2] < 0 then Φ[i-1/2] :Φ[i] ;

patankar notation: (1980)

p10:F[i-1/2]*Φ[i-1/2] = Φ[i-1]*max(F[i-1/2],0) -Φ[i]*max(-F[i-1/2],0);

p10a:F[w]*Φ[w] =Φ[W]*max(F[i-1/2],0) -Φ[p]*max(-F[i-1/2],0);/*see figure 1*/

------------------------------------------------------------------------------------------------------------------

pantakar notation : < F[.],0 > = max(F[.],0), operator : max .
-------------------------------------------------------------------------------------------------------------------