kill(all)$

maxima :convection part no source problem,max function,version 1.0
peter.vlasschaert@gmail.com,10/11/2022

rem : non uniform grid .

load :packages , option
------------------------------

load(lrats)$ ratprint:false$

convection (with max ,upwind scheme)
----------------------------------------------------

ref : see website
---------------------

grid spacing  : different
=================
same:
new:31/10/2022:maxima :  appendix : fvm : Grid generation for CV in 1D,version 1.0
integration : CV
===========
new:20/10/2022:maxima :  fvm: code d(c)integrate(expr,var,lo,hi),version 1.0
upwind        : 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

rem : method introduced by Patankar & Spalding

x-dirction : graph

legend :
------------------------------------------------------------------
||.|| = distance
||pe|| = ||pw|| = Δx
||a1a2|| = ||Wp|| =Δx(w)
||a2a3|| = ||Ep|| =  Δx(e)
=====================================
rem :||aiaj|| = ||ajai|| ,{i,j} = integer (≠ 0)
------------------------------------------------------------------
nodal points : {W,p,E}
interface : {w,e}
-------------------------------------------------------------------
orange : CV = Control Volume ,CS = Control Surface
→ in : w  (CS)
→ out : e (CS)
CV : w ≤ x ≤ e
--------------------------------------------------------------------
pos : direction :
in   :F(w) = u(w)*ρ(w) → Φ (w)
out :F(e)  = u(e)*ρ(e)   → Φ (e)
---------------------------------------------------------------------
interface :{Φ (w),Φ (w)}
nodal      :{Φ (W),Φ (p),Φ (E)}:know
----------------------------------------------------------------------

max function :
===========
reason implimentation exit Fortran 

example 1: working
--------------------------

s1:[max(2,5),max(-2,5),max(2,-5),max(-2,-5)];

example 2: with assume.
--------------------------------

assume(s2>0);

s1a:max(0,-s2);

reason : max function used by Spalding ,Patankar → direction of flow 
x-direction : u ( = velocity), 1D
----------------------------------------
1e) flow u : Left    → Right    (u>0)   : pos-direction
2e) flow u : Right → Left       (u<0)   : neg-direction

problem : advection part ( flow)

1e) 1D model advection (x-direction : u(.),{w,e})
--------------------------------------------------------------

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

2e) integration 1D volume between : w≤ x ≤ e
------------------------------------------------------------

definition : integration

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

dcintegrate(pp2,xx2,aa2,bb2):=block([s11],
s11:integrate(pp2,xx2),
s22:at(s11,xx2=bb2)-at(s11,xx2=aa2),
return(factor(s22))    
)$

p2a:dintegrate(lhs(p1),x,w,e) = rhs(p1);
p2b:dcintegrate(lhs(p1),x,w,e) = rhs(p1);

rem : doesn't work for the convection term.

How to solve problem : "inside differential of p1"
---------------------------------------------------------------

p3a:part(lhs(p1),0);

inside : derivative 
[arguments , variable ,order]

p3b:args(lhs(p1));

→ check : each element derivative or not derivative.
use : p3a, if derivative use  "at"
---------------------------------------------------------------------
rem : dintegrate,dcintegrate
---------------------------------------------------------------------

p3b1:args(p3b[1]);

p3b2:m(p3b[2]);

p3c:args(p3b[1]);

p3c1:makelist(subst([m=p3b1[i]],p3b2),i,1,length(p3b1));

p3c2:product(p3c1[i],i,1,length(p3c1));

integration : [w]→[e]

p3c3:at(p3c2,p3b[2]=e) - at(p3c2,p3b[2]=w);

p3c3a:lratsubst([u(e)*ρ(e) =F(e),u(w)*ρ(w) =F(w)],p3c3);

rem : subst doesn't work.
-------------------------------------------------------------------------------------------------
rem : list constant inside p3c so we can exclude them fröm the integration
→ how : remove them from list p3c
-------------------------------------------------------------------------------------------------

result : replace "d(c)integrate" with more complex version (see website =later).

3e) profile assumption  : upwind scheme .
------------------------------------------------------

Φ[e] = Φ[p] if   F(e) > 0
Φ[e] = Φ[E] if   F(e) < 0

p4a:part(p3c3a,1) =Φ[p]*max(Fe,0) -Φ[E]*max(-Fe,0);

assume(Fe>0);

p4a1:ev(p4a,nouns);

Φ[w] = Φ[W] if   F(w) > 0
Φ[w] = Φ[p]   if   F(w) < 0

p4b:part(p3c3a,2)=-(Φ[W]*max(Fw,0) -Φ[p]*max(-Fw,0));

assume(Fw<0);

p4b1:ev(p4b,nouns);

p4c:p4a1+p4b1;

p4c1:lhs(p4c) =factor(rhs(p4c));

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