Metodi TVD ad alta
risoluzione
Gabriella Puppo
Sommario

Metodi ibridi con flux limiters
 Confronto di diversi limitatori
 Algoritmi di ricostruzione
 Metodi basati su slope limiters
Metodi con flux limiter
Per prevenire le oscillazioni introdotte dagli schemi di ordine
elevato, costruiamo un metodo ibrido conservativo:
Dove F 1 è un flusso numerico del primo ordine (qui useremo il
flusso di Godunov entropico) e F 2 è un flusso numerico del
secondo ordine (qui useremo il flusso di Lax-Wendroff).
Limiter
Il limiter è una funzione Φ che dipende dal misuratore di
regolarità :
Notare che  dipende dalla direzione delle caratteristiche
function hybrid.m
La function hybrid.m implementa un metodo ibrido basato sulla
tecnica del flux limiting.
La function ha:
- un blocco di inizializzazione (come per gli schemi già visti)
- la costruzione dei due flussi numerici F 1 ed F 2
- il calcolo del limiter
- la costruzione del flusso numerico ibrido
- aggiorno la soluzione (come per gli schemi già visti)
- impongo le condizioni al contorno (come per gli schemi già visti)
function [u,x,phi]=hybrid(u0,flux,flux1,t,cfl,limit,ab)
% [u,x]=HYBRID(U0,FLUX,FLUX1,T,CFL,LIMITER) Calcola la soluzione del problema
% u_t+(flux(u))_x=0 su [-1,1] con un metodo ibrido, basato sul
% flusso di Lax Wendroff conservativo e il flusso di
% Godunov entropico
% FLUX e' la funzione di flusso, con derivata FLUX1
% U0 e' il vettore che contiene le condizioni iniziali,T e'
% l'istante finale, CFL e' una stima di max|f'(u)|
% LIMIT e' una variabile stringa, che definisce il limiter
% usato nello schema ibrido.
% I valori possibili sono:
% LIMIT='mm' MinMod limiter
% LIMIT='vl' Van Leer limiter
% LIMIT='cs' constant limiter =0.5
% LIMIT='sb' Super Bee limiter
% [u,x]=HYBRID(U0,FLUX,FLUX1,T,CFL,LIMITER,AB) Calcola la soluzione del
%problema
% u_t+(flux(u))_x=0 sull'intervallo AB=[A,B]
% Le condizioni al contorno sono contenute nella variabile globale
% BC. I valori possibili sono:
% BC='p' Condizioni al contorno periodiche
% BC='f' Condizioni al contorno free-flow
% LANDA0 (global) e' lo scalare t.c. dt=LANDA0*h/CFL
Blocco di inizializzazione
global bc landa0
if nargin < 7
ab=[-1,1];
end
n=length(u0)-2; a=ab(1); b=ab(2);
h=(b-a)/n;
dt=landa0*h/abs(cfl);
nt=fix(t/dt)+1; % arrotonda (T/DT) all'intero immediatamente superiore
dt = t/nt; landa=dt/h;
Calcolo della soluzione:
for kt = 1:nt
% Calcola il flusso di Godunov e quello di LW
f1=fl_godunov(u0,flux,flux1);
[f2,ai]=fl_lw(u0,flux,flux1,landa);
% Calcola il flusso ibrido
phi=limiter(u0,ai,limit);
for i=1:n+1
f(i)=f1(i)+phi(i)*(f2(i)-f1(i));
end
% Calcola la soluzione
for i=2:n+1
u(i) = u0(i) - landa*(f(i)-f(i-1));
end
…. (condizioni al contorno)
end
Ho bisogno di:
- una function che calcoli il flusso di Lax Wendroff
- una function che calcoli il flusso di Godunov
- una function che calcoli il limiter, in modo upwind, cella per cella.
Flusso di Lax-Wendroff: function fl_lw.m
function [f2,ai]=fl_lw(u0,flux,flux1,landa)
% Calcola il flusso numerico con il metodo di Lax Wendroff
% [F2,AI]=FL_LW(U0,FLUX,FLUX1,LANDA)
% F2 = flusso calcolato;
% AI = velocita' caratteristica ai bordi della cella
% U0 = soluzione numerica
% FLUX, FLUX1 = flusso e derivata del flusso
% LANDA = parametro di griglia
n=length(u0)-2;
fl=feval(flux,u0);
for i=1:n+1
uimez(i)=(u0(i)+u0(i+1))/2; % Calcola la media di U
end
ai=feval(flux1,uimez);
% Calcola il flusso di Lax-Wendroff: memorizza f(i+1/2) in f2(i)
for i=1:n+1
f2(i)=(fl(i+1)+fl(i))/2 -landa/2*ai(i)*(fl(i+1)-fl(i));
end
Flusso di Godunov: function fl_godunov.m
function f1=fl_godunov(u0,flux,flux1)
% Calcola il flusso numerico con il metodo di Godunov
% [F1]=FL_GODUNOV(U0,FLUX,FLUX1)
% F1 = flusso calcolato;
% U0 = soluzione numerica
% FLUX, FLUX1 = flusso e derivata del flusso
n=length(u0)-2;
% Flusso di Godunov f1
fl=feval(flux,u0);
fl1=feval(flux1,u0);
for i=1:n+1
s = (fl(i+1) - fl(i))*(u0(i+1)-u0(i));
if s >= 0
f1(i) = fl(i);
else
f1(i) = fl(i+1);
end
anche qui devo aggiungere l’entropy fix
% entropy fix: corregge il flusso se c'e' una rarefazione
% transonica
if fl1(i) <0 & fl1(i+1)> 0
% trova il valore di u per il quale f'(u)=0
if u0(i) >= u0(i+1)
us=fzero(flux1,[u0(i+1),u0(i)]);
else
us=fzero(flux1,[u0(i),u0(i+1)]);
end
f1(i) = feval(flux,us);
end
end
function limiter.m
La function limiter.m ha la seguente struttura:
- Calcola l’indicatore di regolarità in modo upwind
- Calcola il limitatore scegliendo da una lista di limitatori
E’ importante evitare denominatori nulli
Inizializzazione
function phi=limiter(u,ai,stringa)
% PHI=LIMITER(U,STRINGA) calcola il limitatore PHI
% per la funzione di griglia U di tipo STRINGA
% AI = f'(u_(i+1/2)): il segno di AI determina quali
%
punti vengono introdotti nello stencil
% STRINGA = mm (minmod)
%
= cs (costante)
%
= vl (Van Leer)
%
= sb (Super bee)
n=length(u)-2;
Calcolo dell’indicatore di regolarità:
for i=2:n+1
if i==2
du_meno=u(i)-u(i-1);
du_piu=u(i+1)-u(i);
elseif i==n+1
du_meno=u(i)-u(i-1);
du_piu=u(i+1)-u(i);
else
if ai(i) >= 0
du_meno=u(i)-u(i-1);
du_piu=u(i+1)-u(i);
else
du_meno=u(i+1)-u(i);
du_piu=u(i+2)-u(i+1);
end
end
Limitatori Van Leer e MinMod
if stringa=='mm'
if du_meno*du_piu <= 0
phi(i)=0;
elseif abs(du_meno)< abs(du_piu) %& abs(du_piu) > 10.^(-12)
phi(i) = du_meno/du_piu;
else
phi(i) = 1;
end
elseif stringa=='cs'
phi(i)=0.5;
elseif stringa=='vl'
if du_meno*du_piu <= 0
phi(i)=0.0;
else
phi(i)=2*abs(du_meno)/(abs(du_piu)+abs(du_meno));
end
Infine, aggiungo lo script script_hybrid.m che lancia i programmi
% Questo script fa partire il metodo ibrido
% conservativo per equazioni non lineari
% PHI e' il limitatore
clear
global bc landa0
bc='p'; landa0=0.9; n=100;
ab=[-2,2]; t=4; limiter='mm';
u0=init(@quadra,n,ab); cfl=2;
f=inline('x');
f1=inline('1+0*x');
[u,x,phi]=hybrid(u0,f,f1,t,cfl,limiter,ab);
figure
plot(x,u0,'g','Linewidth',2)
hold on
plot(x,u,'b','Linewidth',2)
Risolvo un
problema
di linear advection
con un’onda
quadra come dato
iniziale, con il
limitatore MinMod
Se uso il metodo di Lax Wendroff, ottengo
con le oscillazioni che mi aspetto
Se uso il metodo di Godunov da solo:
che presenta una forte diffusione numerica
Il metodo ibrido con il limitatore MinMod fornisce questo risultato:
che è decisamente meglio.
Con il limitatore di van Leer ottengo:
Notare che la diffusione numerica qui è diminuita
Infine con il limitatore Super Bee ottengo:
Esercizi
Rifare i test precedenti usando una funzione sinusoidale, e
calcolare l’errore in norma 1 e in norma infinito per i diversi
limitatori
Rifare i test precedenti usando l’equazione di Burgers e una
soluzione iniziale di tipo sinusoidale.
Algoritmi di ricostruzione
Un algoritmo di ricostruzione permette di passare dai valori della
soluzione definita per punti su una griglia ad una funzione definita
su tutto un intervallo.
Consideriamo funzioni polinomiali a tratti, usando le medie di
cella della soluzione come dati di interpolazione.
Struttura dei programmi
Per visualizzare il comportamento dell’interpolazione polinomiale
a tratti:
- assegno dei dati U i su una griglia rada, costituita da N intervalli
Ii ;
- calcolo il polinomio di interpolazione su ogni intervallo I i;
- costruisco una griglia fitta su ogni intervallo I i, dove valuto il
polinomio di interpolazione;
- visualizzo i risultati.
Interpolazione costante a tratti
Interpolazione
costante a tratti di
una funzione
regolare.
Notare le
discontinuità ai
bordi delle celle.
Aumentando il numero di nodi della griglia, l’interpolante si
avvicina alla soluzione, e i salti ai bordi delle celle diventano più
piccoli.
Con una griglia ancora più fitta, ottengo
Andamento dell’errore
L’errore diminuisce linearmente, infittendo la griglia:
N=5
N = 10
N = 20
N = 40
err = 0.5878
err = 0.3090
err = 0.1564
err = 0.0785
L’andamento dell’interpolante segue le oscillazioni della soluzione,
senza aumentarle, anche nel caso di una soluzione discontinua:
Interpolazione lineare a tratti
Interpolazione
lineare a tratti
di una
funzione
regolare.
Notare
l’aumento
delle
oscillazioni.
Aumentando il numero di nodi della griglia, l’interpolante si
avvicina rapidamente alla soluzione, e i salti ai bordi delle celle
diventano molto più piccoli.
Andamento dell’errore
L’errore diminuisce quadraticamente, infittendo la griglia:
N=5
N = 10
N = 20
N = 40
err = 0.2163
err = 0.0489
err = 0.0125
err = 0.0031
Inoltre, l’ampiezza delle oscillazioni introdotte dall’interpolante
diminuisce rapidamente, se la soluzione è regolare.
Quando la soluzione presenta delle discontinuità, l’interpolazione
lineare a tratti può essere molto oscillatoria.
Interpolazione quadratica a tratti
Anche in
questo caso,
ottengo dei
buoni risultati
se la soluzione
è regolare
Se invece la soluzione contiene delle singolarità, l’interpolante
quadratico a tratti introduce delle oscillazioni spurie.
Listato dello script per gli algoritmi di
ricostruzione
script_ricostruzione.m
% questo script lancia la routine di ricostruzione
u=inline('sin(pi*x)');
%u=inline('sin(pi*x)+sin(5*pi*x)');
%u=inline('cos(pi*x).*sign(x)');
grado=2;
n=10; ab=[-1,1];
err=ricostruzione(u,grado,n,ab)
Listato della function di ricostruzione
function errore=ricostruzione(u,grado,n,ab)
% ERRORE=RICOSTRUZIONE(U,GRADO,N,AB)
% Disegna il grafico della ricostruzione polinomiale a tratti
% di grado GRADO della funzione U(X) sull'intervallo AB,
% suddiviso in N sottointervalli uguali, e stima l'errore
% ERRORE in norma infinito.
a=ab(1); b=ab(2);
h=(b-a)/n;
x=linspace(a-h/2,b+h/2,n+2);
err=zeros(n,1);
ux=feval(u,x);
figure
% Disegna i punti di interpolazione
plot(x(2:n+1),ux(2:n+1),'r*')
hold on
% Costruisce la griglia su cui calcolare la ricostruzione
dx=h/10;
xx=a;
for j=2:n+1
i=1;
xx=x(j)-h/2;
while x(j)-h/2 -100*eps <= xx & xx <= x(j)+h/2 + 100*eps
if grado==0
% interpola con un polinomio costante a tratti
ug(i)=ux(j);
elseif grado==1
ug(i) = (ux(j+1)-ux(j-1))/(2*h)*(xx-x(j))+ux(j);
elseif grado==2
dd=ux(j+1)-2*ux(j)+ux(j-1);
ug(i) = ux(j) +(ux(j+1)-ux(j-1))/(2*h)*(xx-x(j)) ...
+dd/h^2*(xx-x(j)).^2 -dd/24;
end
xg(i)=xx;
xx=xx+dx; i=i+1;
end
uexa=feval(u,xg);
plot(xg,uexa,'g','Linewidth',2)
plot(xg,ug,'Linewidth',2)
err(j)=norm(ug-uexa,inf);
end
axis([a,b,min(ux)*1.2,max(ux)*1.2])
errore=norm(err,inf);
Esercizi
Calcolare l’andamento dell’errore per l’interpolante
quadratico a tratti per l’interpolazione di una funzione
regolare.
Studiare il comportamento delle oscillazioni spurie presenti
nell’interpolazione di funzioni con gradini, raffinando la
griglia.
Costruire un interpolante lineare a tratti con limitatore di
pendenza
Metodi semidiscreti

Costruzione di un metodo semidiscreto del
secondo ordine
 Calcolo dei valori estrapolati al bordo
 Calcolo dei flussi numerici
 Integrazione nel tempo tramite metodo
Runge-Kutta
Metodo semidiscreto del
secondo ordine
Nel metodo semidiscreto del secondo ordine, risolvo il sistema di
equazioni alle derivate ordinarie con il metodo di Heun:
Per iniziare, identifico le medie di cella ottenute al passo precedente
con il valore iniziale dello schema Runge Kutta
Applico un passo di ricostruzione:
Calcolo il flusso numerico:
Calcolo la soluzione intermedia:
Calcolo i nuovi valori al bordo
Calcolo il nuovo flusso numerico
Infine, calcolo la soluzione al nuovo passo temporale, utilizzando
entrambi i flussi numerici calcolati precedentemente
Quindi devo scrivere:
- una function che calcoli le pendenze, usando un limitatore
- una function che calcoli il flusso numerico utilizzando i due
valori estrapolati al bordo
- una function che implementi il metodo
function slope.m
Questa function calcola le pendenza per la ricostruzione lineare a
tratti
function sigma=slope(u,stringa)
% SIGMA=SLOPE(U,STRINGA) calcola la pendenza SIGMA
% per la funzione di griglia U con il limitatore di tipo STRINGA
% STRINGA = mm (minmod)
%
= cs (costante)
%
= vl (Van Leer)
%
= sb (Super bee)
n=length(u)-2;
for i=2:n+1
du_meno=u(i)-u(i-1);
du_piu=u(i+1)-u(i);
if stringa=='mm'
if du_meno*du_piu <= 0
sigma(i)=0;
elseif abs(du_meno)< abs(du_piu)
sigma(i) = du_meno;
else
sigma(i) = du_piu;
end
elseif stringa=='cs'
sigma(i)=0.5*(u(i+1)-u(i-1));
end
end
% I valori al bordo vengono scelti uguali ai loro vicini
sigma(1)=sigma(2); sigma(n+2)=sigma(n+1);
function fl_god__semidis
function f1=fl_god__semidis(up,um,flux,flux1)
% Calcola il flusso numerico con il metodo di Godunov
% [F1]=FL_GODUNOV(UP,UM,FLUX,FLUX1)
% per i metodi semidiscreti
% F1 = flusso calcolato;
% UP = soluzione numerica in i+1/2 da destra
% UM = soluzione numerica in i+1/2 da sinistra
% FLUX, FLUX1 = flusso e derivata del flusso
% Flusso di Godunov f1
flp=feval(flux,up);
flm=feval(flux,um);
fl1p=feval(flux1,up);
fl1m=feval(flux1,um);
s = (flp - flm).*(up-um);
n=length(up)-2;
for i=1:n+1
if s(i) >= 0
f1(i) = flm(i);
else
f1(i) = flp(i);
end
% entropy fix: corregge il flusso se c'e' una rarefazione
% transonica
if fl1m(i) <0 & fl1p(i)> 0
% trova il valore di u per il quale f'(u)=0
if um(i) <= up(i)
us=fzero(flux1,[um(i),up(i)]);
else
us=fzero(flux1,[up(i),um(i)]);
end
f1(i) = feval(flux,us);
end
end
function semidiscrete.m
function [u,x]=semidiscrete(u0,flux,flux1,t,cfl,limit,ab)
% [u,x]=SEMIDISCRETE(U0,FLUX,FLUX1,T,CFL,LIMITER)
% Calcola la soluzione del problema
% u_t+(flux(u))_x=0 su [-1,1] con un metodo semidiscreto, basato sul
% flusso numerico di Godunov entropico e una ricostruzione
% lineare a tratti con slope limiter.
% Integra nel tempo con il metodo di Heun.
% FLUX e' la funzione di flusso, con derivata FLUX1
% U0 e' il vettore che contiene le condizioni iniziali,T e'
% l'istante finale, CFL e' una stima di max|f'(u)|
% LIMIT e' una variabile stringa, che definisce il limiter
% I valori possibili sono:
% LIMIT='mm' MinMod limiter
% LIMIT='cs' constant limiter =0.5
Blocco di inizializzazione
% Le condizioni al contorno sono contenute nella variabile globale
% BC. I valori possibili sono:
% BC='p' Condizioni al contorno periodiche
% BC='f' Condizioni al contorno free-flow
% LANDA0 (global) e' lo scalare t.c. dt=LANDA0*h/CFL
global bc landa0
if nargin < 7
ab=[-1,1];
end
n=length(u0)-2; a=ab(1); b=ab(2);
h=(b-a)/n;
dt=landa0*h/abs(cfl);
nt=fix(t/dt)+1; % arrotonda (T/DT) all'intero immediatamente superiore
dt = t/nt; landa=dt/h;
Calcolo del primo flusso di Runge -Kutta
for kt = 1:nt
% Calcola i valori estrapolati al bordo, memorizzando i in i+1/2
sigma=slope(u0,limit);
um=u0+sigma*h/2; %soluz in i+1/2 da sinistra
up(1:n+1)=u0(2:n+2)-sigma(2:n+2)*h/2; %soluz in i+1/2 da destra
up(n+2)=up(n+1); % up e um hanno le stesse dimensioni
% Calcola il flusso di Godunov per i valori al bordo
f1=fl_god_semidis(up,um,flux,flux1);
% Calcola la soluzione intermedia
for i=2:n+1
u(i) = u0(i) -landa*(f1(i)-f1(i-1));
end
… aggiunge le condizioni al bordo …
Applica la ricostruzione alla soluzione intermedia e calcola il
secondo flusso di Runge-Kutta
% Calcola i valori estrapolati al bordo
sigma=slope(u,limit);
um=u0+sigma*h/2; %soluz in i+1/2 da sinistra
up(1:n+1)=u0(2:n+2)-sigma(2:n+2)*h/2; %soluz in i+1/2 da destra
up(n+2)=up(n+1); % up e um hanno le stesse dimensioni
% Calcola il flusso di Godunov per i valori al bordo
f2=fl_god_semidis(up,um,flux,flux1);
Infine assembla il flusso globale e aggiorna la soluzione
% Assembla il flusso globale
ftot=0.5*(f1 + f2);
% Calcola la soluzione al passo n+1
for i=2:n+1
u(i) = u0(i) -landa*(ftot(i)-ftot(i-1));
end
… condizioni al contorno…
u0=u; % aggiorna u0 per il prossimo passo
end
x=linspace(a-h/2,b+h/2,n+2);
Script script_semidis.m
Per lanciare il metodo semidiscreto, utilizzo questo script:
Inizializzazione:
% Questo script fa partire il metodo semidiscreto
% conservativo per equazioni non lineari
% PHI e' il limitatore
clear
global bc landa0
bc='p'; landa0=0.9;
Per risolvere un problema di convezione lineare, con un dato
iniziale ad onda quadra
ab=[-2,2]; t=4; limiter='mm'; n=200;
u0=init(@quadra,n,ab); cfl=1;
f=inline('x');
f1=inline('1+0*x');
[u,x]=semidiscrete(u0,f,f1,t,cfl,limiter,ab);
figure
plot(x,u0,'g','Linewidth',2)
hold on
plot(x,u,'b','Linewidth',2)
title('Metodo ibrido')
Per risolvere l’equazione di Burgers con dato iniziale di tipo
inline
ux=inline('1+0.5*sin(pi*x)'); cfl=1.5;
ab=[-2,2]; t=4; limiter='mm';
f=inline('0.5*x.^2');
f1=inline('x');
[u,x]=semidiscrete(u0,f,f1,t,cfl,limiter,ab);
figure
plot(x,u0,'g','Linewidth',2)
hold on
plot(x,u,'b','Linewidth',2)
Linear advection, metodo semidiscreto con MinMod, T=4
Eq. di Burgers, metodo semidiscreto con MinMod, T=0.5
Anche l’onda transonica viene colta bene
Eq. di Burgers, metodo flux limiter con MinMod, T=0.5
Il metodo con flux limiter non e’ entropico.
Esercizi
Considerare un problema di linear advection, con dato iniziale
regolare. Calcolare l’errore usando il metodo flux limiter e il
metodo semidiscreto, con gli stessi parametri computazionali.
Quale metodo da’ l’errore piu’ piccolo?
Considerare un problema di linear advection, con dato iniziale a
gradino. Calcolare la soluzione usando il metodo flux limiter e il
metodo semidiscreto, con gli stessi parametri computazionali.
Quale metodo ha una risoluzione migliore?