1 / 63

Metodi TVD ad alta risoluzione

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.

conner
Download Presentation

Metodi TVD ad alta risoluzione

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Metodi TVD ad alta risoluzione Gabriella Puppo

  2. Sommario • Metodi ibridi con flux limiters • Confronto di diversi limitatori • Algoritmi di ricostruzione • Metodi basati su slope limiters

  3. 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).

  4. Limiter Il limiter è una funzione Φ che dipende dal misuratore di regolarità : Notare che  dipende dalla direzione delle caratteristiche

  5. 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)

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

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

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

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

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

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

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

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

  14. 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;

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

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

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

  18. Se uso il metodo di Lax Wendroff, ottengo con le oscillazioni che mi aspetto

  19. Se uso il metodo di Godunov da solo: che presenta una forte diffusione numerica

  20. Il metodo ibrido con il limitatore MinMod fornisce questo risultato: che è decisamente meglio.

  21. Con il limitatore di van Leer ottengo: Notare che la diffusione numerica qui è diminuita

  22. Infine con il limitatore Super Bee ottengo:

  23. 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 inizialedi tipo sinusoidale.

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

  25. 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 • 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.

  26. Interpolazione costante a tratti Interpolazione costante a tratti di una funzione regolare. Notare le discontinuità ai bordi delle celle.

  27. Aumentando il numero di nodi della griglia, l’interpolante si avvicina alla soluzione, e i salti ai bordi delle celle diventano più piccoli.

  28. Con una griglia ancora più fitta, ottengo

  29. Andamento 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

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

  31. Interpolazione lineare a tratti Interpolazione lineare a tratti di una funzione regolare. Notare l’aumento delle oscillazioni.

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

  33. Andamento dell’errore L’errore diminuisce quadraticamente, infittendo la griglia: N = 5 err = 0.2163 N = 10 err = 0.0489 N = 20 err = 0.0125 N = 40 err = 0.0031 Inoltre, l’ampiezza delle oscillazioni introdotte dall’interpolante diminuisce rapidamente, se la soluzione è regolare.

  34. Quando la soluzione presenta delle discontinuità, l’interpolazione lineare a tratti può essere molto oscillatoria.

  35. Interpolazione quadratica a tratti Anche in questo caso, ottengo dei buoni risultati se la soluzione è regolare

  36. Se invece la soluzione contiene delle singolarità, l’interpolante quadratico a tratti introduce delle oscillazioni spurie.

  37. 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)

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

  39. 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;

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

  41. 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);

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

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

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

  45. Applico un passo di ricostruzione: Calcolo il flusso numerico: Calcolo la soluzione intermedia: Calcolo i nuovi valori al bordo Calcolo il nuovo flusso numerico

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

  47. 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)

  48. 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);

  49. functionfl_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

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

More Related