Metodi TVD ad alta Metodi TVD ad alta risoluzionerisoluzione

Gabriella Puppo

SommarioSommario

Metodi ibridi con flux limitersConfronto di diversi limitatoriAlgoritmi di ricostruzioneMetodi basati su slope limiters

Metodi con flux limiterMetodi con flux limiterPer 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).

LimiterLimiterIl limiter è una funzione Φ che dipende dal misuratore di regolarità :

Notare che dipende dalla direzione delle caratteristiche

function hybrid.mfunction 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

global bc landa0if nargin < 7 ab=[-1,1];endn=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 superioredt = t/nt; landa=dt/h;

Blocco di inizializzazione

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

Calcolo della soluzione:

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.

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 griglian=length(u0)-2;fl=feval(flux,u0);for i=1:n+1 uimez(i)=(u0(i)+u0(i+1))/2; % Calcola la media di Uendai=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 Lax-Wendroff: function fl_lw.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 flusson=length(u0)-2;% Flusso di Godunov f1fl=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

Flusso di Godunov: function fl_godunov.m

% 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); endend

anche qui devo aggiungere l’entropy fix

function limiter.mfunction 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

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;

Inizializzazione

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

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

Limitatori Van Leer e MinMod

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 limitatoreclearglobal bc landa0bc='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);figureplot(x,u0,'g','Linewidth',2)hold onplot(x,u,'b','Linewidth',2)

Risolvo un problemadi 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:

EserciziEsercizi

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 ricostruzioneAlgoritmi 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 programmiStruttura dei programmi

Per visualizzare il comportamento dell’interpolazione polinomiale a tratti:- assegno dei dati U i su una griglia rada, costituita da N intervalli I i ;- 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 trattiInterpolazione 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’erroreAndamento dell’errore

L’errore diminuisce linearmente, infittendo la griglia:

N = 5 err = 0.5878 N = 10 err = 0.3090 N = 20 err = 0.1564 N = 40 err = 0.0785

L’andamento dell’interpolante segue le oscillazioni della soluzione, senza aumentarle, anche nel caso di una soluzione discontinua:

Interpolazione lineare a trattiInterpolazione 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’erroreAndamento dell’errore

N = 5 err = 0.2163N = 10 err = 0.0489N = 20 err = 0.0125N = 40 err = 0.0031

L’errore diminuisce quadraticamente, infittendo la griglia:

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 trattiInterpolazione 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 Listato dello script per gli algoritmi di ricostruzionericostruzione

script_ricostruzione.m

% questo script lancia la routine di ricostruzioneu=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 ricostruzioneListato 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 interpolazioneplot(x(2:n+1),ux(2:n+1),'r*')hold on% Costruisce la griglia su cui calcolare la ricostruzionedx=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);endaxis([a,b,min(ux)*1.2,max(ux)*1.2])errore=norm(err,inf);

EserciziEserciziCalcolare 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 semidiscretiMetodi semidiscreti

Costruzione di un metodo semidiscreto del secondo ordine

Calcolo dei valori estrapolati al bordoCalcolo dei flussi numericiIntegrazione nel tempo tramite metodo

Runge-Kutta

Metodo semidiscreto del Metodo semidiscreto del secondo ordinesecondo 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 precedentecon 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 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)

Questa function calcola le pendenza per la ricostruzione lineare a tratti

functionfunction slope.m slope.m

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)); endend% I valori al bordo vengono scelti uguali ai loro vicinisigma(1)=sigma(2); sigma(n+2)=sigma(n+1);

functionfunction fl_god__semidisfl_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 [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

functionfunction semidiscrete.m semidiscrete.m

% 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 landa0if nargin < 7 ab=[-1,1];endn=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 superioredt = t/nt; landa=dt/h;

Blocco di inizializzazione

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

Calcolo del primo flusso di Runge -Kutta

… aggiunge le condizioni al bordo …

% 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);

Applica la ricostruzione alla soluzione intermedia e calcola il secondo flusso di Runge-Kutta

% 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 endx=linspace(a-h/2,b+h/2,n+2);

Infine assembla il flusso globale e aggiorna la soluzione

ScriptScript script_semidis.m script_semidis.m

% Questo script fa partire il metodo semidiscreto% conservativo per equazioni non lineari% PHI e' il limitatoreclearglobal bc landa0bc='p'; landa0=0.9;

Per lanciare il metodo semidiscreto, utilizzo questo script:

Inizializzazione:

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);figureplot(x,u0,'g','Linewidth',2)hold onplot(x,u,'b','Linewidth',2)title('Metodo ibrido')

Per risolvere un problema di convezione lineare, con un dato iniziale ad onda quadra

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);figureplot(x,u0,'g','Linewidth',2)hold onplot(x,u,'b','Linewidth',2)

Per risolvere l’equazione di Burgers con dato iniziale di tipo inline

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.

EserciziEsercizi

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?