kill(all)$

maxima :diffusion+source problem,a[p],a[E]a[W],b,version 1.0
peter.vlasschaert@gmail.com,10/11/2022

rem : non uniform grid .

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

load(lrats)$ ratprint:false$

convection - diffusion problem ( integration),exact
-------------------------------------------------------------------

grid spacing  : different
=================
same:
new:31/10/2022:maxima :  appendix : fvm : Grid generation for CV in 1D,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
----------------------------------------------------------------------

depends([Tp,T],x);

steady case : heat conduction with source (S = linear equations)
-----------------------------------------------------------------------------------

p1: 'diff(k*'diff(T,x),x)+S=0;

Φ = function ( T = temperature)
Γ = k  (heat )

p1a:[S=Sc+Sp*Tp];

p2:subst(p1a,p1);

1e) integration : between interface ([w]&[e])
------------------------------------------------------------

How to integrate : use  "dcintegrate"

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

a) need integration : variable from differential.

p3a1:args(part(lhs(p2)));

p3a2:args(p3a1[1]);

p3a3:p3a2[2];

b) integration :interface [w] → [e]

p3b:dintegrate((p2),p3a3,w,e);

p3b1:lratsubst([e-w=Δx],p3b);

p3b2:(lhs(p3b1)) = rhs(p3b1);

p3b3:expand(p3b2);

2e) discritization : central  difference.
------------------------------------------------

p4a:p3a2[1];

p4b:makelist(at(p4a*k(i),p3a3=i),i,[w,e])/k;

→ central difference

p4b1:part(p4b[1],1)*k =(T[p]-T[W])/Δx(w)*k(w);
p4b2:part(p4b[2],1)*k =(T[E]-T[p])/Δx(e)*k(e);
p4b3:Tp=T[p];

p4b4:[p4b1,p4b2,p4b3];

p4c:subst(p4b4,p3b3);

p4c1:expand(p4c)*(-1);

p4c1a:a[p] =expand(ratcoeff(lhs(p4c1),T[p]));

p4c1b:T[p]*a[p] =T[p]*rhs(p4c1a);

p4c1c:expand(lhs(p4c1)-rhs(p4c1b));

other side : equation 

p4c1d:p4c1c*(-1);

3e) final result
--------------------

see website :
=========
---------------------------------------------------------------------------------------------------------------
new:09/11/2022:maxima :  convection -diffusion problem (exact : solution),version 1.0
---------------------------------------------------------------------------------------------------------------

p5a:a[p]*T[p] =a[E]*T[E] +a[W]*T[W]+b;

p5a1:rhs(p5a) = p4c1d;

p5a2:[b,T[W],T[E]];

p5a3:makelist(ratcoeff(p5a1,p5a2[i]),i,[2,3]);

p5aa3:makelist(reverse(p5a3[i]),i,1,length(p5a3));

p5b3:[subst(makelist(p5a2[i]=0,i,[2,3]),p5a1)];

p5c3:[lratsubst(p5aa3,p4c1a)];

p5d3:flatten(makelist([p5a3,p5b3,p5c3]));

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