kill(all)$

maxima : fvm :implementation (TVD) for 1d convection-diffusion problem,version 1.0 
peter.vlasschaert@gmail.com,26/11/2022

load packages
------------------

load(lrats)$

implementation (TVD)
=================

see website :
-----------------
ref 1:
new:20/10/2022:maxima :  fvm: code d(c)integrate(expr,var,lo,hi),version 1.0
download  (maxima,code)        :maximafvmdintegrateverison12022
ref 2:
new:21/11/2022:maxima : :TVD framework,find r, Ψ(r) ,  version 1.0
download  (maxima,code)       :maximafvmTVDframeworkschemeversion12022final1
ref 3:
maxima : fvm :convection -diffusion relation peclet cell number (max :limit),version 1.0
download  (maxima,code)       :maximaconvdiffusionandrelationtopecletcelnumnversion12022final

rem : CDS = Central Difference Scheme (for diffusion)

fvm:finite volume method.

convection - diffusion equation 
-----------------------------------------

s1:'diff(u(x)*ρ(x)*Φ(x),x) = 'diff(Γ(x)*diff(Φ(x),x),x);

integration expr contains : differential  , between : aa1 ≤ xx1 ≤ bb1:
--------------------------------------------------------------------------------------

dintegrate(pp1,xx1,aa1,bb1):=block([s1],
s1:integrate(pp1,xx1),
s2:at(s1,xx1=bb1)-at(s1,xx1=aa1),
return(s2)    
)$

s1a:dintegrate(lhs(s1),x,w,e) =dintegrate(rhs(s1),x,w,e);

discretised version s1a
-----------------------------

1e) convection term : TVD scheme (lhs) ,ref 2
2e) diffusion     term : CDS scheme (rhs)
3e) combine 1e) & 2e)

→ convection term:

step 1:

p2a1:flatten(makelist([F[v]=u(v)*ρ(v),Φ[v]=Φ(v)],v,[e,w]));

p2a2:makelist(rhs(p2a1[v])=lhs(p2a1[v]),v,1,length(p2a1));

p2a3:lratsubst(p2a2,lhs(s1a));

step 2: "TVD ,x-direction  ( u > 0)"

find : Φ[f] , faces (e,w) , finite volume [w & e]
------------------------------------------------------------
index U = upwind       → Φ[U]
          D = downwind  → Φ[D]
use: more General discription takes in account the "direction velocity of u ".

ref2: {p2b4,p5b}

p2a4a:Φ[f]=(Φ[D]-Φ[U])·Ψ(r)/2+Φ[U];

p2a4b:Ψ(r)=(Φ[UU]-Φ[U])/(Φ[U]-Φ[D]);

part : u > 0
----------------

nodal values : 1D spatial 
-------------------------------
u > 0 ,   UU > U > p > D > DD  (= WW>W>w> p >e > E > EE), because u > 0
u < 0,    DD > D > p > U > UU  (= WW>W>w> p >e> E  > EE)
-------------------------------------------------------------------------
E   : east  -> EE
W : west -> WW
-------------------------------------------------------------------------
general : U,D

p2a4c:subst([p2a4b],p2a4a);

step 3: face {w,e}

→ face:w ,see : nodal values. ( w = ref),see also Figure 1 "above"

p2a4bw:[f=w,r=w,U=W,UU=WW,D=p];

p2a4bw1:subst(p2a4bw,p2a4b);

p2a4bw2:subst(p2a4bw,p2a4a);

p2a4bwf:subst([p2a4bw1],p2a4bw2);

→ face:e ,see : nodal values. ( e = ref),see also Figure 1 "above"

p2a4be:[f=e,r=e,U=p,UU=W,D=E];

p2a4be1:subst(p2a4be,p2a4b);

p2a4be2:subst(p2a4be,p2a4a);

p2a4bef:subst([p2a4be1],p2a4be2);

step 4: input into (p2a3)

p2a51:subst([p2a4bw2,p2a4be2],p2a3);

p2a52:subst([p2a4bw1,p2a4be1],p2a51);

→ diffusion term:  CDS scheme (ref 3)

use ref 3: p3a

p2b1:[at('diff(Φ(x),x,1),x=w)=(Φ[p]-Φ[W])/Δx(w),at('diff(Φ(x),x,1),x=e)=(Φ[E]-Φ[p])/Δx(e)];

rhs : s1a

p2b2:rhs(s1a);

use :p2b1

p2b3a:subst(p2b1,p2b2);

set : D[f] = Γ(f)/Δx(f) , f ={w,e}

p2b3b1:makelist(Γ(f)/Δx(f)=D[f],f,[w,e]);

p2b3b2:factorsum(lratsubst(p2b3b1,p2b3a));

rem : factor , collectterms
https://maxima.sourceforge.io/docs/manual/maxima_324.html#collectterms
-------------------------------------------------------------------------------------------------------------------------------------------------------
Function: factorfacsum (expr, arg_1, ... arg_n)
Returns a form of expr which is obtained by calling facsum on the factors of expr with arg_1, ... arg_n as arguments.   
If any of the factors of expr is raised to a power, both the factor and the exponent will be processed in this way. 
-------------------------------------------------------------------------------------------------------------------------------------------------------

→ combine

use :{p2a51,p2b3b2}

p2c1:p2a51=p2b3b2;

general equation :

p3:a[p] *Φ[p] = a[W]*Φ[W] + a[E]*Φ[E]+S[TVD];

find : S[TVD]

p3a1:subst([Ψ(e)=0,Ψ(w)=0],lhs(p2c1));

p3a2:S[TVD]=collectterms(expand(lhs(p2c1)-p3a1),F[w],Ψ(w),F[e],Ψ(e));

find :{a[p],a[W],a[E]}

use : p3a1

p3a3:p3a1=rhs(p2c1);

→ a[W]

p3a3b1:ratcoeff(p3a3,Φ[W]);

p3a31a:a[W]=rhs(p3a3b1)-lhs(p3a3b1);

→ a[E]

p3a3b2:a[E]=ratcoeff(rhs(p3a3),Φ[E]);

→ a[p]

p3a3b3:ratcoeff(p3a3,Φ[p]);

p3a3b3:a[p]=lhs(p3a3b3)-rhs(p3a3b3);/*reverse : rhs*/

use :a[W],a[E]

p3a3b3a:lratsubst([reverse(p3a31a),reverse(p3a3b2)],p3a3b3);

part : u < 0
----------------

nodal values : 1D spatial 
-------------------------------
u > 0 ,   UU > U > p > D > DD  (= WW>W>w> p >e > E > EE), 
u < 0,    DD > D > p > U > UU  (= WW>W>w> p >e> E  > EE), because u < 0
-------------------------------------------------------------------------
E   : east  -> EE
W : west -> WW
-------------------------------------------------------------------------
general : U,D

→ 1e) step  (TVD)

face {w,e}

→ face [e], e = ref

p4a1e:[f=e,r=e,U=E,UU=EE,D=p];

p4a1e1:subst(p4a1e,p2a4b);

p4a1e2:subst(p4a1e,p2a4a);

p4a1e3:subst([p4a1e1],p4a1e2);

→ face [w], w = ref

p4a1w:[f=w,r=w,U=p,UU=E,D=W];

p4a1w1:subst(p4a1w,p2a4b);

p4a1w2:subst(p4a1w,p2a4a);

p4a1w3:subst([p4a1w1],p4a1w2);

→ 2e) step : "convection -diffusion " equation

use :p2a3

p5a1:p2a3;

p5a2:subst([p4a1e2,p4a1w2],lhs(p5a1))=p2b3b2;

→ 3e) step : {S[TVD],a[p],a[W],a[E]}

S[TVD]

p5a3:subst([Ψ(e)=0,Ψ(w)=0],lhs(p5a2));

p5a3b:p5a3= rhs(p5a2);

p5a31:S[TVD] =collectterms(expand(lhs(p5a2)-p5a3),F[w],Ψ(w),F[e],Ψ(e));

general equation :

p5bb:a[p] *Φ[p] = a[W]*Φ[W] + a[E]*Φ[E]+S[TVD];

=>  a[E],use : p5a3b

p5be:ratcoeff(p5a3b,Φ[E]);

p5be1:a[E]=rhs(p5be)-lhs(p5be);

=>  a[W],use : p5a3b

p5bw:a[W]=rhs(ratcoeff(p5a3b,Φ[W]));

=>  a[p],use : p5a3b

p5bpp:ratcoeff(p5a3b,Φ[p]);

p5bp1:a[p]=lhs(p5bpp)-rhs(p5bpp);/*reverse : rhs*/

p5bp2:lratsubst([reverse(p5be1),reverse(p5bw)],p5bp1);

final result
---------------

→ TVD neighbourhood
---------------------------------

p7a:a[W] =D[w]+max(F[e],0);
p7b:a[E] =D[e]+max(-F[w],0);

→ TVD(c),correction source term.
---------------------------------------------
α(:) =1 → F(.) > 0   {. = e,w}
α(:) =0 → F(.) <0

p7c:S[TVD(α)]=1/2*F[e]*((1-α(e))*Ψ(e)-α(e)*Ψ(e))*(Φ[E]-Φ[p])+1/2*F[w]*((1-α(w))*Ψ(w)-α(w)*Ψ(w))*(Φ[p]-Φ[W]);

p7:a[p] =rhs(p7a)+rhs(p7b)+(F[e]-F[w]);

use : conservation law. in → out (F[e] =F[w]),no internal production (CV)
-------------------------------------------------------------------------------------------------
example : chemical reaction,source inside cv.

p7cl:a[p] =rhs(p7a)+rhs(p7b);

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