kill(all)$

maxima:fvm: TVD framework,find r, Ψ(r) ,  version 1.0 
peter.vlasschaert@gmail.com,21/11/2022

TVD Framework
-----------------------

TVD = Total Variation Diminshing

rem : TV for convective discretized solution 
i = index refer to nodal values : WW,W,p,E,EE (=spatial domain)

definition : TV
-------------------

p1:TV = sum(abs(Φ[i+1]-Φ[i]),i,1,n);

n = nodal values in the spatial domain 

Definition : TVD
----------------------
numerical scheme TVD =TV solution dimishes in the "time"

rem : solution be stable and non-oscillatory scheme preserving : monotonicity
( = create non- local extrema)
→ doesn`t create overshoot and undershoot
-----------------------------------------------------------------------------------------------------
result : TV must not increase

p1a:TV(Φ^(m+1))  ≤  TV(Φ^(m));

rem : m & m+1 are successive time steps.

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".

1e : step UDS = Upwind Difference Scheme
---------------------------------------------------------

new:04/11/2022:maxima :  appendix : fvm : upwind-scheme Φ[w],version 1.0 when(u>0,u<0),version 1.0
new:04/11/2022:maxima :  appendix : fvm : upwind-scheme Φ[e],version 1.0 when(u>0,u<0),version 1.0

p2a:Φ[f] =Φ[U];

2e : step CDS = Upwind Difference Scheme (  = mean value scheme )
---------------------------------------------------------

p2b:Φ[f] =1/2*(Φ[U]+Φ[D]);

p2b1:lhs(p2b) = rhs(p2a) +1/2*nn;

→ find nn ?

p2b2:solve(rhs(p2b)=rhs(p2b1),nn);

p2b2a:p2b2[1];

→ final result:

p2b3:subst(p2b2,p2b1);

p2b3a:rhs(p2b3);

p2b3b:ratcoeff( p2b3a,rhs(p2b2a))*(-1);

use :Ψ(r)

p2b4:lhs(p2b)=rhs(p2a)+Ψ(r)*p2b3b*rhs(p2b2a);

step 3: find Ψ(r) , UDS & CDS
---------------------------------------

→ UDS ,use p2a

p3c1:rhs(p2a) = rhs(p2b4);

p3c1a:solve(p3c1,Ψ(r))[1];

→ CDS ,use p2b3

p3c2:rhs(p2b) = rhs(p2b4);

p3c2a:solve(p3c2,Ψ(r))[1];

result :Φ[U]=((Φ[D]-Φ[U])*Ψ(r))/2+Φ[U]
----------------------------------------------------
UDS : Ψ(r)=0
CDS : Ψ(r)=1
----------------------------------------------------

step 4: QUICK scheme
------------------------------

see : ref

new:04/11/2022:maxima :  appendix : fvm : quickscheme Φ[w],version 1.0 when(u>0,u<0),version 1.0
new:04/11/2022:maxima :  appendix : fvm : quickscheme Φ[e],version 1.0 when(u>0,u<0),version 1.0

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

quick scheme : U,UU,D 
------------------------------

p4:Φ[f] =Φ[U]+1/8*(3*Φ[D]-2*Φ[U]-Φ[UU]);

see ref :

p4a:ratsimp(expand(p4));

check :8 = 6+3-1

use : p3c1

p4b1: p3c1;

find :Ψ(r)

p4b2:rhs(p4) =rhs(p4b1);

p4b2a:p4b2-Φ[U];

p4b2b:factor(solve(p4b2a,Ψ(r))[1]);

set : r =(Φ[U]-Φ[UU])/(Φ[D]-Φ[U])

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

inside : p4b2b

p4c2:lratsubst([reverse(p4c1)],p4b2b);

result :Φ[f]=((Φ[D]-Φ[U])*Ψ(r))/2+Φ[U]
----------------------------------------------------
QUICK : Ψ(r)=(r+3)/4
----------------------------------------------------
set r=(Φ[U]-Φ[UU])/(Φ[D]-Φ[U])

step 5: other scheme 
------------------------------

→  SOU scheme relation with Φ[f]=((Φ[D]-Φ[U])*Ψ(r))/2+Φ[U]

p5:Φ[f] =Φ[U]+1/2*(Φ[U]-Φ[UU]);/*SOU scheme*/
p5a:p2b4;
p5b:solve(rhs(p5)=rhs(p2b4),Ψ(r))[1];
p5c:lratsubst(reverse(p4c1),p5b);

→FROMM scheme

p5d:Ψ(r) = (1+r)/2;

prove :Sweby [1984]
---------------------------
→ scheme TVD (necessary and sufficient condition)
1e)=>  0<r<1 → upperlimit Ψ(r) =2*r → TVD : Ψ(r) ≤ 2r
2e)=>      r≥1 → upperlimit Ψ(r) =2    → TVD : Ψ(r) ≤ 2

rem : above are all linear schemes.
----------------------------------------------
(UDS,CDS,SOU,FROMM,QUICK)

result : Convective Schemes
=====================

→ linear schemes
-----------------------------------------------------------
UDS         : Ψ(r) = 0 (very stable highly difusive)
CDS         : Ψ(r) = 1 (highly dispersive               )
SOU         : Ψ(r) = r
FROMM   : Ψ(r) = (1+r)/2
QUICK     : Ψ(r) = (3+r)/4
-------------------------------------------------------------
→ point intersection (r & Ψ(r))
r=1,Ψ(r) =1 → (CDS,SOU,FROMM,QUICK)
--------------------------------------------------------------

→ non-linear schemes
-----------------------------------------------------------
Min Mod : Ψ(r) = max(0,min(1,r))
Van Leer : Ψ(r) = (r+abs(r))/(1+abs(r))
MUSCL   : Ψ(r) = max(0,min(2*r,0.5*r+0.5,2))
SMART   : Ψ(r) = max(0,min(2*r,0.75*r+0.25,4))
-------------------------------------------------------------

step 6: plot scheme 
-----------------------------

wxplot2d([2,2*r], [r,0,5],[y,0,4],
 [nticks,5])$

blue line :Ψ(r) = 2
red   line :Ψ(r) = 2*r

TVD regio : Monotonicity
---------------------------------

red-line to intersection(blue line,red line) and below the blue-line

H(r) := if r < 1 then 2*r else 2 $

wxplot2d ([H(r)], [r, 0, 5], [y, 0, 2.5],[xlabel,"r"],[ylabel,"Ψ(r)"],[title,"TVD"])$

below : blue-line

ref : use plot2d
===========
-----------------------------------------------------------------------
(%i1) xy: [[10, .6], [20, .9], [30, 1.1], [40, 1.3], [50, 1.4]]$
(%i2) plot2d([[discrete, xy], 2*%pi*sqrt(l/980)], [l,0,50],
        [style, points, lines], [color, red, blue],
        [point_type, asterisk],
        [legend, "experiment", "theory"],
        [xlabel, "pendulum's length (cm)"],
        [ylabel, "period (s)"])$
-----------------------------------------------------------------------

r-Ψ(r),wxplot2d (= Different Convective Scheme )
-----------------------------------------------------------------
→ UDS,CDS,SOU,QUICK
→ intersection : "see above:result : Convective Schemes"

wxplot2d ([H(r),0,1,r,(3+r)/4], [r, 0, 5.2], [y, -0.1, 2.2],[xlabel,"r"],[ylabel,"Ψ(r)"],[title,"TVD"])$

output :
----------

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

UDS    :       TVD
CDS    : not TVD  : r <0.5
SOU    : not TVD  : r > 2
QUICK: not TVD  : r<3/7 & r>5 

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

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