POLITECNICO DI MILANOCorso di Laurea in Ingegneria Matematica

Codice 2D ai volumi finiti per unproblema di Darcy

Relatore: Prof. Luca Formaggia

Tesi di Laurea di:An-phi Nguyen, matricola 757314

Anno Accademico 2012-2013

Ai miei genitori

Indice

1 Introduzione 11.1 Obiettivo e struttura della tesi . . . . . . . . . . . . . . . . . . 1

2 Il problema di Darcy 22.1 Modello fisico-matematico . . . . . . . . . . . . . . . . . . . . 2

2.1.1 Le leggi di conversazione . . . . . . . . . . . . . . . . . 22.1.2 La legge di Darcy . . . . . . . . . . . . . . . . . . . . . 4

2.2 Modello numerico . . . . . . . . . . . . . . . . . . . . . . . . . 52.2.1 Ipotesi semplificative . . . . . . . . . . . . . . . . . . . 5

3 Stato dell’arte 63.1 Metodo alle differenze finite . . . . . . . . . . . . . . . . . . . 63.2 Metodo agli elementi finiti . . . . . . . . . . . . . . . . . . . . 73.3 Metodi a volumi finiti . . . . . . . . . . . . . . . . . . . . . . 9

4 Il Control Volume Finite Element Method 10

5 Algoritmo in C++ 155.1 Cenno alle librerie e ai tool utilizzati . . . . . . . . . . . . . . 155.2 File .vtk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165.3 Classi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185.4 Il metodo load() . . . . . . . . . . . . . . . . . . . . . . . . . . 23

5.4.1 Cancellazione dei dati precedenti . . . . . . . . . . . . 235.4.2 Lettura dell’intestazione . . . . . . . . . . . . . . . . . 235.4.3 Lettura dei punti, delle celle e dei lati . . . . . . . . . . 245.4.4 Lettura degli attributi e calcolo dei flussi numerici . . . 255.4.5 Costruzione della matrice e del vettore dei termini noti 25

6 Esempi di utilizzo 286.1 Un esempio . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

6.1.1 Installazione della libreria . . . . . . . . . . . . . . . . 28

3

6.1.2 Descrizione dell’esempio . . . . . . . . . . . . . . . . . 296.1.3 Esecuzione e analisi dei risultati . . . . . . . . . . . . . 296.1.4 Visualizzazione grafica dei risultati . . . . . . . . . . . 31

6.2 Altri esempi . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316.2.1 Esempio a dominio generico . . . . . . . . . . . . . . . 326.2.2 Esempio con permeabilita non costante . . . . . . . . . 34

7 Conclusione 36

A Analisi dei coefficienti di trasmissibilita 38A.1 Principio del massimo e condizioni sufficienti . . . . . . . . . . 38A.2 Interpretazione geometrica . . . . . . . . . . . . . . . . . . . . 40

B Eccezioni 41

C La classe Point2D 43

D Esempio 44

4

Capitolo 1

Introduzione

La Fluidodinamica e la parte della meccanica relativa alla dinamica dei liq-

uidi e dei gas1. Essa si occupa di costruire modelli di fenomeni naturali odi applicazioni industriali coinvolgenti fluidi. Troppo spesso, pero, trovareuna soluzione analitica del problema risulta molto difficile o, addirittura, im-possibile. D’altra parte, compiere esperimenti potrebbe risultare oltremododispendioso. Per questo motivo, negli ultimi decenni, si e sviluppata la Flu-idodinamica Computazionale (CFD) che si occupa della risoluzione numeri-ca, con l’ausilio dei computer, dei suddetti modelli. Un esempio di problemastudiato nell’ambito di questa disciplina e la simulazione di flussi in mezziporosi.

1.1 Obiettivo e struttura della tesi

Il presente lavoro di tesi ha come scopo quello di presentare e fornire unasoluzione numerica per un modello di flussi di Darcy in mezzi porosi. Perfare cio, e stata sviluppata una libreria C++ che verra presentata piu indettaglio nei prossimi capitoli.L’esposizione sara suddivisa principalmente in tre parti:

1. Presentazione del problema

2. Stato dell’arte e presentazione del metodo risolutivo

3. Struttura dell’algoritmo

Infine verranno presentati brevi esempi di utilizzo (con relative istruzioni)2.

1da Enciclopedia Treccani2Come sistema operativo sara utilizzato Linux.

1

Capitolo 2

Il problema di Darcy

I flussi in mezzi porosi3 sono un modello fluidodinamico di grande importan-za in molte applicazioni, quali lo studio delle falde acquifere o delle riservepetrolifere. I motivi per cui questo tipo di analisi e cosı importante si possonotrovare sulla maggior parte dei libri riguardante questo argomento, come peresempio [Fanchi]. In questo ambito, notevole fu il contributo di Henry Darcy(ingegnere francese, 1803-1858), da cui prende il nome il modello matematicopiu utilizzato per tale tipo di flussi. Argomento del seguente capitolo saran-no le nozioni fondamentali alla base dell’algoritmo implementato. Per unatrattazione piu dettagliata si vedano ad esempio [Pinder] o [Chavent].

2.1 Modello fisico-matematico

Per la costruzione di un modello matematico occorrono due componenti fon-damentali: leggi generali e relazioni costitutive. In questo caso, le leggi gen-erali sono le tipiche leggi di conservazione e di bilancio della Meccanica deiContinui e la relazione costitutiva e la legge (sperimentale) di Darcy. Persemplificare la trattazione, si rendono necessarie alcune ipotesi: i fenomenidi dispersione possono essere trascurati e la superficie di contatto tra fluidoe solido e impermeabile al fluido, cioe non puo essere attraversata. Inoltre,verranno considerati solo fluidi monofase.

2.1.1 Le leggi di conversazione

Siano x = (x1, x2, x3) e t le variabili spaziali e temporali rispettivamente.Siano, inoltre, φ la porosita (cioe la frazione di volume occupabile dal fluido)

3Per mezzo poroso si intende un materiale contenente dei pori che permettono iltransitare di fluidi attraverso il materiale stesso.

2

del mezzo, ρ la densita del fluido per unita di volume, u = (u1, u2, u3) la ve-locita di Darcy4 e q le sorgenti (o pozzi) esogene. Si consideri un volumetto

Figura 2.1: Volume differenziale

(Fig 2.1) con spigoli paralleli agli assi coordinati e centro (x1, x2, x3). Glispigoli nella direzione dell’i-esimo asse hanno lunghezza ∆xi. Il flusso dimassa per unita di superficie nella direzione i e dato da ρui. Ora, facendoriferimento alla figura 2.1, la massa totale entrante nel volume infinitesi-mo nella direzione 1 e (ρu1)x1−

∆x12

∆x2∆x3, mentre la massa totale uscente e

(ρu1)x1+∆x12

∆x2∆x3. Analogamente si ottengono i flussi di massa nelle altre

direzioni.La velocita di variazione totale di massa nel volumetto e data da ∂(φρ)

∂t∆x1∆x2∆x3.

L’accumulo (o la “distruzione”) di massa ad opera di sorgenti (o pozzi) lo-calizzati nel volume e data da q∆x1∆x2∆x3.Facendo un bilancio di massa, otteniamo la seguente espressione:

(ρu1)x1−∆x12

∆x2∆x3 − (ρu1)x1+∆x12

∆x2∆x3+

(ρu2)x2−∆x22

∆x1∆x3 − (ρu2)x2+∆x22

∆x1∆x3+

(ρu3)x3−∆x32

∆x1∆x2 − (ρu1)x3+∆x32

∆x1∆x2

=∂(φρ)

∂t∆x1∆x2∆x3 − q∆x1∆x2∆x3

4La velocita di Darcy e la velocita media di un fluido attraverso a una sezione del mezzotrasversale al flusso.

3

Ora, dividendo per ∆x1∆x2∆x3 e passando al limite per ∆xi → 0 ∀i, siottiene5:

∂(φρ)

∂t+∇ · (ρu) = q

2.1.2 La legge di Darcy

In ogni punto x ∈ R3 e definito un vettore u, detto velocita di Darcy, tale

che la portata volumetrica del fluido attraverso una superficie ds di normaleν e data da u · νds.La legge di Darcy ha la seguente forma6

u = −K(x)

µ[∇p− ρg∇z(x)]

dove:

• K(x) e la permeabilita al fluido in questione del mezzo poroso nel puntox. Per mezzi isotropi K e uno scalare, mentre per mezzi anisotropi e untensore. In particolare, si puo dimostrare che e un tensore simmetricodefinito positivo. La permeabilita e misurata in Darcy [L]2.

• µ e la viscosita dinamica del fluido, misurata in Poise [M ][L]−1[s]−1.

• p e la pressione del fluido in Pascal [M ][L]−1[s]−2.

• ρ e la densita [M ][L]−3.

• g e l’accelerazione di gravita [L][s]2.

• z(x) e la quota del punto in esame [L].

5∇· e l’operatore di divergenza. Dato un campo vettoriale v ∈ C1

∇ · v =∂v1

∂x1

+∂v2

∂x2

+∂v3

∂x3

6∇ e l’operatore gradiente. Dato un campo scalare T ∈ C1:

∇T = (∂T

∂x,∂T

∂y,∂T

∂z)

4

2.2 Modello numerico

Riassumendo i risultati del precedente paragrafo, le equazioni che caratter-izzano il problema in esame sono:

∂(φρ)∂t

+∇ · (ρu) = q

u = −K(x)µ

[∇p− ρg∇z(x)]

2.2.1 Ipotesi semplificative

Nella seguente trattazione ci limiteremo a un caso particolare del problemadi Darcy. Le ipotesi che si utilizzeranno sono:

• stazionarieta

• fluido omogeneo (ρ costante)

• variazioni di quota trascurabili rispetto alle variazioni di pressione

Di conseguenza, il problema da studiare e (inglobando il coefficiente di vis-cosita in quello di permeabilita):

∇ · u = f

u = −K(x)∇pInoltre verra considerato solo il caso bidimensionale.

5

Capitolo 3

Stato dell’arte

Spesso accade, soprattutto in ambito matematico, che i problemi possanoessere risolti in molti modi e quello presentato in questo lavoro di tesi non fasicuramente eccezione. Di seguito verranno presentati alcuni metodi comune-mente usati nella risoluzione di modelli per flussi di Darcy7. Come problemadi riferimento si consideri quello esposto al capitolo 2.

3.1 Metodo alle differenze finite

L’idea alla base di questo approccio e quello di approssimare le derivatespaziali e temporale mediante i rispettivi rapporti incrementali.Comunemente una derivata (parziale) prima viene approssimata in uno deiseguenti 3 modi:

1. In avanti

∂p

∂x1

(x1, x2) =p(x1 + h, x2)− p(x1, x2)

h

2. All’indietro

∂p

∂x1

(x1, x2) =p(x1, x2)− p(x1 − h, x2)

h

3. Centrato

∂p

∂x1

(x1, x2) =p(x1 + h, x2)− p(x1 − h, x2)

2h

7Per una descrizione piu approfondita si veda [Chen]

6

dove h e, in questo caso, l’incremento spaziale nella direzione x1.Si consideri ora un dominio di calcolo Ω = [a; b] × [c; d] rettangolare e lo sisuddivida in rettangoli piu piccoli. L’obiettivo e quello di calcolare le pres-sioni nei vertici di questi rettangolini8. I nodi cosı individuati sono identificatidalle coordinate (x1,i, x2,j) con 0 ≤ i ≤ n e 0 ≤ j ≤ m

x1,i ∈ [a; b] e x2,j ∈ [c; d]

x1,0 = a, x1,n = b, x2,0 = c, x2,m = d

x1,i− 1

2=

1

2(x1,i−1 + x1,i)

h1,i = x1,i − x1,i−1

Tenendo a mente le discretizzazioni della derivata prima presentate a iniziosessione, si puo mostrare che l’equazione alle derivate parziali

−∇ · (K(x)∇p) = f

si puo approssimare con

−

Ki+ 1

2,jpi+1,j−pi,j

h1,i+1

2

−Ki− 1

2,jpi,j−pi−1,j

h1,i− 1

2

h1,i

−

Ki,j+ 1

2

pi,j+1−pi,jh2,j+1

2

−Ki,j− 1

2

pi,j−pi,j−1

h2,j− 1

2

h2,j

= fi,j

dove i pedici indicano i punti della griglia in cui e valutata la permeabilita9.Cio porta a dover risolvere un sistema di equazioni della forma

−Ki+ 1

2,j(pi+1,j−pi,j)+Ki− 1

2,j(pi,j−pi−1,j)−Ki,j+ 1

2(pi,j+1−pi,j)+Ki,j− 1

2(pi,j−pi,j−1) = Fi,j

dove Fi,j = fi,jh1,ih2,j.Le condizioni al bordo si impongono modificando opportunamente la matricedei coefficienti e il termine noto del sistema lineare ottenuto10.

3.2 Metodo agli elementi finiti

Questo metodo utilizza un tipo di discretizzazione diverso dal precedente. Percomprenderlo meglio, si consideri il problema semplificato monodimensionale:

8Questo particolare tipo di griglia e detto point-distributed grid. In alternativa si potevaconsiderare una griglia in cui i nodi fossero rappresentati dai centri dei suddetti rettangoli.

9Come visto precedentemente, K(x) e in generale un tensore del second’ordine. Pernon appesantire la notazione, e stato preso in considerazione solo il caso di K(x) scalare.

10Per esempio, le condizioni di Dirichlet si possono imporre portando a secondo membrole pressioni note, cioe quelle appartenenti alla frontiera del dominio.

7

−p′′(x) = f(x)p(0) = p(1) = 0

con 0 < x < 1.Si suddivida l’intervallo [0; 1] in M+1 intervalli, cioe si scelgano M+2 puntitali che 0 < x0 < x1 < x2 < ... < xM < xM+1 = 1. Si definiscono diconseguenza gli incrementi hi = xi − xi−1 con i = 1...M + 1. Si definisca lospazio lineare finito

Vh = v : v e una funzione continua su [0; 1] ,dv

dxe

costante a tratti e v(0) = v(1) = 0

dotato del prodotto scalare

(v, w) =

∫ 1

0

v(x)w(x)dx

Si puo dimostrare che le funzioni ϕi ∈ Vh con i = 1...M tali che ϕi(xj) = δijcostituiscono una base di Vh. Inoltre ogni v ∈ Vh ha una rappresentazioneunica della forma

∑M

i=1 v(xi)ϕi(x). Di conseguenza, moltiplicando entrambi imembri dell’equazione a derivate parziali di partenza per una generica v ∈ Vh,integrando per parti (tenendo conto delle condizioni iniziali) e sostituendo ap la sua approssimazione ph(x) =

∑M

i=1 p(xi)ϕi(x) ∈ Vh, si ottiene

M∑

i=1

pi(dϕi

dx,dv

dx) = (f, v)

con pi = p(xi). Ma poiche ogni v ∈ Vh puo essere espresso come combinazionelineare delle funzioni costituenti la base, ci si puo limitare a risolvere il sistemalineare

M∑

i=1

aijpi = fj

con j = 1, ...,M , aij = (dϕi

dx,dϕj

dx) e fj = (f, ϕj).

Tutto il ragionamento esposto puo essere esteso al caso bidimensionale ot-tenendo ancora un sistema lineare. Si puo inoltre dimostrare che la matricedei coefficienti ottenuta e in entrambi i casi simmetrica e definita positiva11.

11Una matrice n× n simmetrica A si dice definita positiva se xTAx > 0 ∀x ∈ Rn

8

3.3 Metodi a volumi finiti

Questi metodi possono essere considerati una sorta di variante dei metodiagli elementi finiti.Infatti, essi mantengono la stessa idea di approssimare la soluzione tramitefunzioni arbitrarie che, nell’approccio piu semplice, sono combinazioni di poli-nomi di primo grado. La differenza sta nel fatto che questa discretizzazioneviene effettuata su un volume di controllo costruito intorno ai nodi della meshdi calcolo. Un esempio e l’algoritmo che prende il nome di Control Volume

Finite Element Method (CVFE) ed e l’algoritmo implementato nel presentelavoro di tesi. Pertanto si rimandano i dettagli del metodo al prossimo capi-tolo per una trattazione piu approfondita.Il metodo CVFE puo essere generalizzato utilizzando, ad esempio, funzionispline12. Questa generalizzazione prende il nome di Control Volume Finite

Approximation Method.

12Dati dei punti (detti nodi) appartenenti al dominio della funzione da interpolare, perfunzione spline si intende una funzione costituita da un insieme di polinomi raccordati neinodi in modo tale da essere, in questi nodi, continua e con derivata continua fino a unordine prestabilito

9

Capitolo 4

Il Control Volume FiniteElement Method

Come gia accennato nel capitolo precedente, il metodo Control Volume Finite

Element e una variante del metodo a elementi finiti. Si consideri un dominio

Figura 4.1: Volume di controllo

limitato Ω ⊂ R2 e il problema

−∇ · (K(x)∇p) = f in Ω

L’algoritmo che si intende implementare gestisce solo mesh di calcolo a ele-menti triangolari. Dato un dominio qualsiasi bisogna quindi ottenerne unatriangolazione. Questo aspetto del problema esula dallo scopo della presentetrattazione. Per approfondimenti si consulti [Quarteroni].Siano Mh la suddetta mesh e M i relativi elementi. Sia Vi un volume dicontrollo. I Vi sono costruiti intorno ai nodi della mesh unendo i punti medi

10

dei lati dei triangoli con un punto qualsiasi al loro interno (Fig. 4.1). Sos-tituendo la pressione p con ph ∈ Vh = v ∈ C(M) : v lineare su ∀M ∈ Mhe integrando su Vi si ottiene:

−

∫

Vi

∇ · (K(x)∇ph)dx =

∫

Vi

fdx

Per il teorema della divergenza:

−

∫

∂Vi

K(x)∇ph · νdl =

∫

Vi

fdx (4.1)

Concetriamoci ora sul termine a sinistra dell’uguale.Si puo osservare che ∇ph · ν e continuo attraverso ogni segmento di ∂Vi (chegiace all’interno di un triangolo). Di conseguenza, se K(x) e continuo, losara anche il flusso K(x)∇ph · ν. Inoltre, (4.1) indica che questo metodo elocalmente conservativo.Dato un triangolo M con vertici mi,mj ,mk, punti medi ma,mb,md e centromc (Fig 4.2). L’approssimazione ph di p su M e data da

ph = p(mi)λi + p(mj)λj + p(mk)λk (4.2)

dove λi sono le funzioni nodali, tali che λi(mj) = δij con

Figura 4.2: Un triangolo M della mesh

λi + λj + λk = 1 (4.3)

Si definiscanoai = mj,2 −mk,2

bi = −(mj,1 −mk,1)

ci = mj,1mk,2 −mj,2mk,1

11

dove mi = (mi,1,mi,2) e i, j, k sono permutati ciclicamente. Si puo provareche le funzioni nodali sono date da

λi

λj

λk

=1

2|M |

ci ai bicj aj bjck ak bk

1x1

x2

(4.4)

dove |M | e l’area del triangolo M. Di conseguenza,

∂λl

∂x1

=al

2|M |,∂λl

∂x2

=bl

2|M |(4.5)

con l = i, j, k.Calcoliamo, ora, (4.1) sul tratto mamcmd (Fig 4.2).Lungo mamc risulta

ν =((mc,2 −ma,2), (ma,1 −mc,1))

|mamc|,

mentre lungo mcmd

ν =((md,2 −mc,2), (mc,1 −md,1))

|mcmd|.

Da cui, tenendo conto di (4.2) e (4.5), si ottiene

−

∫

mamc∪mcmd

K(x)∇ph · νds = |M |k

∑

l=i

K∇λl · ∇λipl (4.6)

Considerando (4.3) e (4.6), si giunge alla seguente formulazione

−

∫

mamc∪mcmd

K(x)∇ph · νds = −Tij(pj − pi)− Tik(pk − pi) (4.7)

dove Tij e Tik sono i coefficienti di trasmissibilita e sono pari a

Tij = −|M |K∇λj · ∇λi

Tik = −|M |K∇λk · ∇λi

Considerando che ogni segmento che collega due nodi (di cui almeno unointerno alla mesh) e adiacente a due triangoli M1 e M2, il coefficiente totaledi trasmissibilita tra due nodi i e j e

Tij = −2

∑

l=1

(|M |K∇λj · ∇λi)

∣

∣

∣

∣

Kl

(4.8)

12

Percio l’equazione lineare13 per il volume di controllo Vi centrato sul nodo i

e

−∑

j∈Ωi

Tij(pj − pi) =

∫

Vi

fdx (4.9)

dove Ωi e l’insieme dei nodi collegati al nodo i.Passiamo ora ad analizzare il termine di destra di (4.1).Sotto l’ipotesi che i termini di sorgente siano costanti sugli elementi dellamesh, risulta

∫

M

fdx = f

∫

M

dx = f |M | (4.10)

con M elemento della mesh Mh.Si osservi che per la costruzione del volume di controllo Vi era stato sceltoun punto qualsiasi interno agli elementi M nell’intorno del nodo i. Si prendaquindi il baricentro in modo tale che ogni M contribuisca alla misura |Vi| perun terzo della sua misura |M |. Grazie a questa scelta, il termine di destra in(4.1) e (4.9) si puo riscrivere

∫

Vi

fdx =∑

k∈Mh,i

fk

3|Mk| (4.11)

Da cui, la (4.9) puo essere finalmente scritta nella sua forma finale

−∑

j∈Ωi

Tij(pj − pi) =∑

k∈Mh,i

fk

3|Mk| (4.12)

con Mh,i ⊂ Mh insieme degli elementi della mesh nell’intorno del nodo i.Iterando il ragionamento su ogni nodo, si costruisce il sistema lineare la cuirisoluzione fornira i valori delle pressioni nei nodi.Pero occorre prima imporre le condizioni al contorno. In questa trattazionesono state considerate solo due tipi di condizioni.

• Condizioni di Dirichlet

p(x) = g(x) ∀x ∈ ∂Ω

Nella risoluzione numerica, questo vincolo si traduce nella versionediscreta

pi = p(xi) = g(xi) ∀i tale che xi ∈ ∂Ω

13Per avere risultati soddisfacenti dal punto di vista fisico, i coefficienti di trasmissi-bilita devono essere non-negativi, in caso contrario la soluzione potrebbe discostarsi dallasituazione reale. Per un approfondimento di questa condizione si veda l’Appendice A.

13

Il modo piu semplice di implementare questa condizione e quella diridurre i gradi di liberta del sistema legati ai nodi noti grazie ai datial bordo. Operativamente, i pi noti devono essere portati a secondomembro nella (4.9) e considerati termini noti del sistema lineare.

• Condizioni di Neumann omogenee

∇p(x) · n(x) = 0 ∀x ∈ ∂Ω

dove n(x) indica un versore normale a ∂Ω. In realta questa condizione equella implementata naturalmente con la costruzione del sistema (4.9)e pertanto non servono altri accorgimenti per imporre questo tipo didati al bordo.

14

Capitolo 5

Algoritmo in C++

In questo capitolo saranno esposte le parti fondamentali dell’algoritmo le cuibasi teoriche sono state presentate nel Capitolo 1.

5.1 Cenno alle librerie e ai tool utilizzati

Per la creazione del programma sono state utilizzate alcune librerie gia es-istenti.

• Standard Library

La STL e stata inclusa per poter adoperare alcune sue funzionalitamolto utili per l’implementazione dell’algoritmo. Ne sono un esempiole eccezioni, gli smart pointers e i container.

• Boost (Versione 1.0)E stato fatto uso della libreria Boost per quanto riguarda la conversionedi stringhe in formati numerici.

• Eigen (Versione 3.1.2)Si e usufruito della libreria Eigen per la gestione dei tensori di perme-abilita e della matrice globale di trasmissibilita e per risolvere il sistemalineare ad essa associato.

Inoltre, sono stati adoperati altri tool:

• CMake

Questo tool e servito per la creazione dei makefile atti all’installazionedella libreria e per la compilazione dei file di esempio14.

14Si veda il Capitolo 6.

15

• Matlab

L’ambiente Matlab e stato utilizzato per implementare script atti allagenerazione di alcuni file esempio, in quanto fornisce un valido algorit-mo di triangolazione di domini.

• Paraview

Valido software di visualizzazione. Con esso si possono visionare siai dati in ingresso al programma che la soluzione fornita. Per unadescrizione del formato di entrambi i file si veda il prossimo paragrafo.

5.2 File .vtk

Per una maggiore comprensione dei capitoli successivi, verra brevemente es-posta la struttura dei file forniti in input/output al/dal programma.Il formato dei file in questione e vtk 15, un formato largamente utilizzato inambito scientifico per la rappresentazione dei dati.Il file e strutturato in tre parti.

• Intestazione

In questa sezione sono fornite le informazioni per decodificare corret-tamente il file. Essa puo essere ulteriormente suddivisa in 3 parti.

1 # vtk DataFi le Vers ion 1 .02 2D Unstructured Grid o f Linear Tr i ang l e s3 ASCII

La prima riga riporta la versione del DataFile.La seconda parte consta di 256 caratteri terminati con un caratterenullo ’\0’. In questo spazio e possibile inserire una descrizione deldataset. La terza parte indica il formato del file: ASCII o BINARY.NOTA:l’algoritmo gestisce solo file in ASCII.

• Struttura dei dati

Questa sezione descrive la geometria dei dati.

1 DATASET UNSTRUCTUREDGRID2

3 POINTS n type4 p0 x p0 y p0 z5 . .6 p(n−1) x p(n−1) y p(n−1) z

15Per la documentazione ufficiale si consulti http://www.vtk.org/VTK/help/documentation.html

16

7

8 CELLS n c e l l n i n f9 n l a t i p0 1 p0 2 . . p0 ( n l a t i )

10 . . .11

12 CELL TYPES n c e l l13 type014 type115 . . .

Innanzitutto viene esplicitato il tipo di dataset. In questo caso si escelto di adoperare una griglia non strutturata16 in quanto per il fun-zionamento del codice non si richiede che i nodi seguano un particolareordine.Successivamente, viene fornito il dataset vero e proprio.Per i punti (cioe i nodi) bisogna specificarne il numero e il tipo con cuivengono descritti (es. float). E da notare che bisogna riportare tutte etre le coordinate dei punti: poiche l’algoritmo tratta solo mesh 2D leterze componenti saranno poste pari a zero.Per ogni cella (cioe gli elementi finiti della mesh), la prima informazioneda fornire e il numero di lati. Gli interi successivi (sulla linea) riportanol’indice dei punti che costituiscono i vertici della cella (ordinati in sen-so antiorario). Il numero totale di celle e di interi per descriverle tuttesono riportati (rispettivamente come n cell e n inf ) nell’intestazionedi questa sottosezione.Nell’ultima parte sono riportati i tipi di celle17. Per il funzionamentodi questo algoritmo, tutti i valori dovranno essere pari a 5.

• Attributi del dataset

Questa sezione fornisce informazioni fondamentali per la costruzionedella matrice di trasmissibilita.

1 POINT DATA n po int2 SCALARS incogn i t a i n t3 LOOKUPTABLE de f au l t4 . .5

6 SCALARS pr e s s i on e f l o a t7 LOOKUPTABLE de f au l t8 . .

16Per una descrizione piu dettagliata dei tipi di dataset utilizzabili si faccia riferimentoalla documentazione ufficiale.

17Si faccia riferimento alla documentazione ufficiale.

17

9

10 CELL DATA n c e l l11 SCALARS sorgente f l o a t12 . .13

14 TENSORS permeab i l i t y f l o a t15 . .

Per quanto riguarda i nodi, i primi dati da riportare sono degli int (chevanno intesi come booleani) per indicare se il valore di pressione neipunti e incognito. Successivamente, bisogna fornire i valori di pressioneper i nodi noti. Per gli altri possono essere inseriti dei valori qualsiasiin quanto non verranno letti. Infine per le celle, occorre fornire lasorgente/pozzo e i tensori di permeabilita.

Il file generato in output dal programma non differisce dal file input. Laprincipale differenza sta nel fatto che, nella sezione relativa agli attributi deldataset, vengono riportati i valori di pressione calcolati dall’algoritmo,per inodi, e i vettori velocita di Darcy, per gli elementi finiti della mesh.

1 POINT DATA n po int2 SCALARS pr e s s i on e f l o a t3 LOOKUPTABLE de f au l t4 . .5

6 CELL DATA n c e l l7 VECTORS v e l o c i t y f l o a t8 . .

5.3 Classi

La logica dell’algoritmo si basa soprattutto sulle formule (1.8) e (1.12). In-fatti, da esse si puo constatare quali siano le informazioni importanti daconsiderare nella costruzione della matrice globale (caratteristiche geomet-riche della mesh, come i vertici dei triangoli, i gradienti delle funzioni nodalietc..., e quelle fisiche del mezzo poroso studiato). Alla luce di cio, sono stateimplementate 4 nuove classi18.

18Le dichiarazioni di queste classi sono contenute nell’header MeshCVFE.hpp e sonoracchiuse nel namespace MeshCvfe. Per poter far uso della libreria e necessaria l’opzione(al momento della compilazione) [$ -lMeshCVFE]

18

• NodoCVFE

1 typede f index type i n t ;2

3 c l a s s NodoCVFE: pub l i c Point2D4 pr i va t e :5 double p r e s s i on e =0;6 bool i n cogn i t a=true ;7 index type id ;8 index type id unk ;9 pub l i c :

10 NodoCVFE( double X, double Y, double p , index type index ) :Point2D (X,Y) , p r e s s i on e (p) , i n cogn i t a ( f a l s e ) , id (

index ) 11 NodoCVFE( double X, double Y, index type index ) : Point2D (

X,Y) , id ( index ) 12 index type showIndex ( ) const r e turn id ;13 bool unknown ( ) const r e turn in cogn i t a ;14 double p r e s su r e ( ) const r e turn p r e s s i on e ;15 index type showUnkIndex ( ) const r e turn id unk ;16 void a s s i gnPr e s su r e ( double p) p r e s s i on e=p ; i n cogn i t a=

f a l s e ; r e turn ;17 void setUnkId ( index type id ) id unk=id ; r e turn ;

:18 f r i e nd bool operator <( const NodoCVFE& a , const NodoCVFE

& b) r e turn a . id<b . id ;19 f r i e nd bool operator==(const NodoCVFE& a , const

NodoCVFE& b) r e turn a . id==b . id ;20 f r i e nd ostream& operator<<(ostream&, const NodoCVFE&) ;21 ;

Questa classe implementa l’unita base della mesh: i nodi. La primacaratteristica da osservare e che essa e una sottoclasse di Point2D.Quest’ultima rappresenta un punto in R

2 e ne contiene le coordinate(x, y).NodoCVFE aggiunge ad essa informazioni importanti per il metodo,quali la pressione, un valore booleano che indica se la pressione nelnodo e incognita o e nota, un indice globale e un indice di incogni-ta19 per semplificare l’interazione con il vettore soluzione del sistemalineare legato al metodo. Un aspetto interessante da considerare e ladefinizione del tipo dell’indice. Attualmente l’indice e stato implemen-tato come int. Solitamente il range del tipo int e sufficiente per indiciz-zare tutti i nodi di una mesh. Naturalmente puo essere reimplementato

19La numerazione di entrambi gli indici parte da 0.

19

con altri tipi quali long int o unsigned int. La scelta e ricaduta su int

poiche risulta essere il piu efficiente, anche se in realta, per l’uso che sifa di esso in questo programma, la scelta del tipo di intero non impattasulla performance in modo significativo. Altri punti significativi sonogli operatori comparativi (l’ordinamento si basa sugli indici globali deinodi), la cui utilita sara chiara piu avanti.

• ElementoCVFE

1 typede f shared ptr<const NodoCVFE> cptrNodo ;2

3 c l a s s ElementoCVFE4 pr i va t e :5 vector<cptrNodo> v e r t i c i ;6 double area=0;7 double so rgente =0;8 MatrixXd Kp;9 VectorXd i ;

10 VectorXd j ;11 VectorXd k ;12 void ca l co l aArea ( ) ;13 pub l i c :14 ElementoCVFE( cptrNodo nodo i , cptrNodo nodo j , cptrNodo

nodo k ) ;15 ElementoCVFE( cptrNodo nodo i , cptrNodo nodo j , cptrNodo

nodo k , Matrix2d&& K) ;16 void s e tPe rmeab i l i t y ( const MatrixXd& K) Kp=K; return ;17 void s e tPe rmeab i l i t y (Matrix2d&& K) Kp=move(K) ; r e turn ;18 const MatrixXd& showPermeabi l i ty ( ) const r e turn Kp;19 const VectorXd& showVectorGrad ( const NodoCVFE& nodo )

const ;20 const vector<cptrNodo>& showVertex ( ) const r e turn

v e r t i c i ;21 double showArea ( ) const r e turn area ;22 ;

La classe ElementoCVFE rappresenta gli elementi finiti della mesh.Contiene un vettore a nodi costanti (vertici), in quanto deve poteravere accesso alle informazioni relative ad essi ma senza il rischio dimodificarli. I vertici sono ordinati in senso anti-orario. Grazie ad essi,viene calcolata l’area tramite la formula dell’area di Gauss20. Kp e iltensore di permeabilita considerato costante su tutto l’elemento. I tre

20Formula di Gauss: A =∑

n

i=1

∣

∣

∣

∣

xiyi+1−yixi+1

∣

∣

∣

∣

con n numero di vertici e la convenzione

n+ 1 = 1.

20

vettori i, j, k sono i vettori gradiente della funzioni nodali dell’elemen-to.Si noti che i vettori gradiente e la matrice di permabilita sono statiimplementati come vettori e matrici dinamiche in previsione di unpossibile ampliamento del codice al caso tridimensionale.

• EdgeCVFE

1 typede f shared ptr<const NodoCVFE> cptrNodo ;2 typede f shared ptr<const ElementoCVFE> cptrElem ;3 c l a s s EdgeCVFE4 pr i va t e :5 cptrNodo i ;6 cptrNodo j ;7 double f l u s s o i j =0;8 i n t numElAdiacenti=0;9 cptrElem e l 1=nu l lp t r , e l 2=nu l l p t r ;

10 pub l i c :11 EdgeCVFE( cptrNodo a , cptrNodo b) ;12 void addElem( cptrElem ) ;13 void c a l c o l a ( ) ;14 double showFlusso ( ) const r e turn f l u s s o i j ;15 cptrNodo showNodoI ( ) const r e turn i ;16 cptrNodo showNodoJ ( ) const r e turn j ;17 f r i e nd bool operator <( const EdgeCVFE& a , const EdgeCVFE

& b) r e turn ( (∗ ( a . i )==∗(b . i ) ) ? (∗ ( a . j )<∗(b . j ) ) : ( ∗ ( a. i )<∗(b . i ) ) ) ;

18 f r i e nd ostream& operator<<(ostream&, const EdgeCVFE&) ;19 ;

Questa classe implementa i lati (o spigoli nel caso 3D) che collegano duenodi della mesh. Caratteristica importante di EdgeCVFE e la gestionedei nodi (a cui si ha accesso tramite puntatore): il nodo i ha indiceglobale < del nodo j. Anche qui, e di particolare interesse l’operatoredi comparazione (che, come si puo notare, fa uso degli operatori < e== della classe nodo). L’utilita di questo operatore e dell’ordinamentodei nodi del lato sara piu chiaro in seguito.I due puntatori a ElementoCVFE tengono conto degli elementi finitiadiacenti al lato.

• MeshCVFE

Questa e la classe che, di fatto, implementa l’algoritmo presentato nelCapitolo 1.

1 typede f shared ptr<ElementoCVFE> ptrElem ;

21

2 typede f shared ptr<EdgeCVFE> ptrEdge ;3 typede f Eigen : : SparseMatrix<double> spar s e ;4

5

6 c l a s s MeshCVFE7 pr i va t e :8 vector<ptrNodo> nodi ;9 vector<ptrElem> e l em e n t i f i n i t i ;

10 vector<ptrEdge> l a t i ;11 s t r i n g MeshINFO ;12 bool so lved=f a l s e ;13 spar s e T;14 VectorXd b ;15 VectorXd so l u z i on e ;16 pub l i c :17 MeshCVFE( ) ;18 MeshCVFE( const s t r i n g f i l ename ) load ( f i l ename ) ;19 index type showNumNodi ( ) r e turn nodi . s i z e ( ) ;20 index type showNumElem( ) r e turn e l em e n t i f i n i t i . s i z e ( )

;21 index type showNumEdge ( ) r e turn l a t i . s i z e ( ) ;22 void s o l v e ( ) ;23 void load ( const s t r i n g f i l ename ) ;24 void printVtk ( const s t r i n g name=” so l u z i on e . vtk” ) const

;25 f r i e nd ostream& operator<<(ostream&,MeshCVFE&) ;26 const NodoCVFE& operator [ ] ( index type i ) const ;27

28 ;

I tre campi fondamentali di questa classe sono: il vettore di nodi, il vet-tore di elementi finiti e il vettore di lati. In realta, non sono contenitoriper le classi vere e proprie, ma per smart pointers. Questo per un piufacile utilizzo condiviso dei dati contenuti in esse, evitando cosı copieinutili delle medesime informazioni. Inoltre, gli smart pointers dan-no un notevole aiuto per quanto riguarda l’allocazione dinamica dellamemoria. Per quanto riguarda gli EdgeCVFE, in fase di lettura e statoscelto di raccoglierli in un set provvisorio per una questione di efficienzanell’inserimento di nuovi lati21. Da qui l’importanza dell’operatore <

della classe EdgeCVFE.La matrice di trasmissibilita e stata implementata come matrice sparsa22.

21Una spiegazione piu dettagliata di cio sara data nella descrizione della funzione load()di questa classe.

22E tipico dei metodi numerici per la soluzione di equazioni a derivate parziali come

22

Altra caratteristica particolare di questa classe e l’implementazionedell’operatore [] che permette un accesso piu rapido ai nodi della mesh.

5.4 Il metodo load()

Di seguito, verra descritto in dettaglio il metodo load() della classe MeshCVFEpoiche costituisce il fulcro dell’intero codice (e la funzione che costruisce ilsistema lineare da risolvere per trovare le pressioni). In ingresso bisognafornire una stringa contenente il nome del file .vtk da aprire (e se necessarioil percorso per arrivare ad esso).Il costruttore della classe MeshCVFE che riceve in ingresso una stringa fun-ziona allo stesso modo (di fatto chiama il metodo load()).Per una presentazione piu chiara, la funzione sara esposta in piu parti.

5.4.1 Cancellazione dei dati precedenti

Il metodo inizia cancellando tutti i dati su cui si e lavorato precedentementee creando un set provvisorio dei lati.

1 // Cance l l a z i one ev en tua l i da t i p r e c eden t i e c r e a z i one de l s e tp r ovv i s o r i o

2 nodi . c l e a r ( ) ;3 e l em e n t i f i n i t i . c l e a r ( ) ;4 l a t i . c l e a r ( ) ;5 so lved=f a l s e ;6 set<ptrEdge , bool (∗ ) ( ptrEdge , ptrEdge )> s l a t i ( compareEdgePointer )

; // s e t p r ovv i s o r i o per i l a t i

5.4.2 Lettura dell’intestazione

Questa sezione del codice si limita a leggere l’intestazione del file .vtk.

1 // I n t e s t a z i o n e2 // Skipping in fo rmat ion about Vers ion3 cout<<”Reading ve r s i on . . . ”<<endl ;4 g e t l i n e ( i n p u t f i l e , t emp st r ing ) ;5 //Get Mesh Informat ion6 cout<<”Reading i n f o . . . ”<<endl ;7 g e t l i n e ( i n p u t f i l e ,MeshINFO) ;8 //Check Format9 cout<<”Reading Data Format . . . ”<<endl ;

quelli presentati nel capitolo 3.

23

10 g e t l i n e ( i n p u t f i l e , t emp st r ing ) ;11 i f ( t emp st r ing==”BINARY” ) throw ErroreASCII ( ) ;12 //Dataset type13 whi le ( i n p u t f i l e >>t emp st r ing ) 14 i f ( t emp st r ing==”DATASET” ) 15 i n p u t f i l e >>t emp st r ing ;16 i f ( t emp st r ing==”UNSTRUCTUREDGRID” ) 17 break ;18 e l s e throw ErroreGrid ( ) ;19 20

Le informazioni riguardanti la versione del file vengono semplicemente igno-rate poiche inutili per il funzionamento del codice. La descrizione del datasetviene salvata nel campo MeshINFO (una string) della mesh. Per il forma-to e la struttura dati viene fatto un controllo per verificare che siano quellirichiesti dal programma. In caso contrario, vengono lanciate delle eccezioni23.

5.4.3 Lettura dei punti, delle celle e dei lati

Da qui si inizia a far largo uso della funzione lexical cast della libreriaboost. Nel caso ci sia qualche errore nella formattazione del file di input, e diconseguenza non venga letto un valore numerico, un’eccezione viene lanciata.Si riporta solo la parte relativa alla lettura delle celle (e dei lati) poiche e lapiu interessante.

1 cout<<”Reading c e l l s . . . ”<<endl ;2 whi le ( i n p u t f i l e >>t emp st r ing ) 3 i f ( t emp st r ing==”CELLS” ) 4 i n p u t f i l e >>t emp st r ing ;5 num els=boost : : l e x i c a l c a s t <double>( t emp st r ing ) ;6 break ;7 8 9 e l em e n t i f i n i t i . r e s e r v e ( num els ) ;

10 i n p u t f i l e >>t emp st r ing ; // sk ipp ing11 f o r ( index type i =0; i<num els ; i++)12 i n p u t f i l e >>t emp st r ing ;13 i n t num lat i ;14 num lat i=boost : : l e x i c a l c a s t <int >( t emp st r ing ) ;15 i f ( num lat i !=3) throw ErroreLat i ( ) ;16 index type temp1 , temp2 , temp3 ;17 i n p u t f i l e >>t emp st r ing ;18 temp1=boost : : l e x i c a l c a s t <index type >( t emp st r ing ) ;

23Per una lista completa delle eccezioni create si veda Appendice B.

24

19 i n p u t f i l e >>t emp st r ing ;20 temp2=boost : : l e x i c a l c a s t <index type >( t emp st r ing ) ;21 i n p u t f i l e >>t emp st r ing ;22 temp3=boost : : l e x i c a l c a s t <index type >( t emp st r ing ) ;23 ptrElem e l (new ElementoCVFE( nodi [ temp1 ] , nodi [ temp2 ] , nodi [

temp3 ] ) ) ;24 e l em e n t i f i n i t i . push back ( e l ) ;25 ptrEdge edge (new EdgeCVFE( nodi [ temp1 ] , nodi [ temp2 ] ) ) ;26 auto r i f=l a t i . i n s e r t ( edge ) ;27 (∗ r i f . f i r s t )−>addElem( e l ) ;28 edge . r e s e t (new EdgeCVFE( nodi [ temp2 ] , nodi [ temp3 ] ) ) ;29 r i f=l a t i . i n s e r t ( edge ) ;30 (∗ r i f . f i r s t )−>addElem( e l ) ;31 edge . r e s e t (new EdgeCVFE( nodi [ temp3 ] , nodi [ temp1 ] ) ) ;32 r i f=l a t i . i n s e r t ( edge ) ;33 (∗ r i f . f i r s t )−>addElem( e l ) ;34 35 l a t i . r e s e r v e ( s l a t i . s i z e ( ) ) ;36 copy ( s l a t i . begin ( ) , s l a t i . end ( ) , i n s e r t e r ( l a t i , l a t i . begin ( ) ) ) ;

Per come e costruito il file di input, e inevitabile che ogni lato (ad eccezionedi quelli appartenenti alla frontiera del dominio di calcolo) sia letto due volte(o piu nel caso di estensione a un caso 3D). Di conseguenza, per evitare diavere dei doppioni, i lati sono implementati temporaneamente come set (perquesto motivo la classe EdgeCVFE e stata dotata di un operatore <). Lafunzione set :: insert(), nel caso che un lato sia gia presente, si limita arestituire il puntatore ad esso.Alla fine di questa fase il set viene convertito a vector<> per velocizzarnel’accesso.

5.4.4 Lettura degli attributi e calcolo dei flussi numeri-ci

In questa sezione (di cui non si riporta il codice poiche non particolarmentesignificativo) vengono letti gli attributi relativi ai punti e alle celle. Comenella fase di lettura dell’intestazione, nel caso di errori nella formattazionedel file, vengono lanciate delle eccezioni.Terminata questa fase, il set dei lati viene attraversato per poter calcolaretutti i flussi numerici (tramite la formula (1.8)).

5.4.5 Costruzione della matrice e del vettore dei ter-mini noti

25

1 cout<<”Bui ld ing Matrix . . . ”<<endl ;2 T. r e s i z e ( unk count , unk count ) ;3 b=std : : move(VectorXd : : Zero ( unk count ) ) ;4 index type i =1;5 index type num edges=l a t i . s i z e ( ) ;6 f o r ( auto i t=l a t i . begin ( ) ; i t != l a t i . end ( ) ; i t++)7 cout<<”Analyzing edge data ”<< i++<<” o f ”<<num edges<<endl ;8 index type i =(∗ i t )−>showNodoI ( )−>showUnkIndex ( ) ;9 index type j =(∗ i t )−>showNodoJ ( )−>showUnkIndex ( ) ;

10 double f =(∗ i t )−>showFlusso ( ) ;11 i f ( (∗ i t )−>showNodoI ( )−>unknown ( ) ) 12 T. c o e f fRe f ( i , i )+=f ;13 double source =0;14 i f ( (∗ i t )−>showNumElAdiacenti ( )==1)15 source=(∗ i t )−>showEl1 ( )−>showSource ( ) ∗(∗ i t )−>showEl1 ( )−>

showArea ( ) /6 ;16 e l s e 17 source =(((∗ i t )−>showEl1 ( )−>showSource ( ) ∗(∗ i t )−>showEl1 ( )−>

showArea ( ) )+((∗ i t )−>showEl2 ( )−>showSource ( ) ∗(∗ i t )−>showEl2 ( )−>showArea ( ) ) ) /6 ;

18 19 b( i )+=source ;20 i f ( (∗ i t )−>showNodoJ ( )−>unknown ( ) ) 21 T. c o e f fRe f ( j , j )+=f ;22 T. i n s e r t ( i , j )=−f ;23 T. i n s e r t ( j , i )=−f ;24 b( j )+=source ;25 e l s e 26 b( i )+=f ∗ ( (∗ i t )−>showNodoJ ( )−>pre s su r e ( ) ) ;27 28 e l s e 29 i f ( (∗ i t )−>showNodoJ ( )−>unknown ( ) ) 30 T. c o e f fRe f ( j , j )+=f ;31 double source =0;32 i f ( (∗ i t )−>showNumElAdiacenti ( )==1)33 source=(∗ i t )−>showEl1 ( )−>showSource ( ) ∗(∗ i t )−>showEl1 ( )−>

showArea ( ) /6 ;34 e l s e 35 source =(((∗ i t )−>showEl1 ( )−>showSource ( ) ∗(∗ i t )−>showEl1 ( )

−>showArea ( ) )+((∗ i t )−>showEl2 ( )−>showSource ( ) ∗(∗ i t )−>showEl2 ( )−>showArea ( ) ) ) /6 ;

36 37 b( j )+=(f ∗ ( (∗ i t )−>showNodoI ( )−>pre s su r e ( ) )+source ) ;38 39 40

26

Queste linee di codice implementano la formula (1.12). Nella matriceglobale di trasmissibilita vengono inseriti i coefficienti relativi alle pressioniincognite, mentre al termine noto contribuiscono le pressioni note e i terminidi sorgente. Si noti che a differenza della formula (1.12), i contributi datidalle sorgenti/pozzi sono divisi per 6 invece che per 3. Questo perche, percome e costruito l’algoritmo, attraversando il set dei lati, ogni triangolo M

della mesh contribuisce due volte ad ogni volume di controllo Vi.

27

Capitolo 6

Esempi di utilizzo

In questa sezione verra mostrato un esempio di funzionamento, insieme a bre-vi istruzioni di installazione e di visualizzazione dei risultati. Infine verrannofatte delle considerazioni sul progetto svolto.

6.1 Un esempio

Verra ora presentato un esempio tratto da [Chen]. Prima forniremo breviistruzioni per l’installazione della libreria. Come detto in precedenza, e nec-essario che siano installati anche le librerie Boost (V 1.0) e Eigen (V 3.1.2).Per maggiori informazioni si consultino le documentazioni ufficiali relative.

6.1.1 Installazione della libreria

Per installare la libreria, si proceda in questo modo:

1. Decomprimere il file CVFE.tar.gz .

2. Se ritenuto opportuno, modificare le variabili LIBRARY OUTPUT

DIRECTORY e HEADER OUTPUT DIRECTORY nel file CMake-Lists.txt per cambiare la directory in cui verranno installati, rispetti-vamente, la libreria e il relativo header file.

3. Creare una cartella in cui allocare i file di output di CMake e accedervi.

4. Lanciare CMake da questa directory [$ cmake <directory in cui sono

stati decompressi i file>]

5. A questo punto sara presente un makefile. Lanciare il commando make.

28

6. Per installare la libreria e l’header, utilizzare il comando make install.Potrebbe essere necessario avere i privilegi di root [$ sudo make install].

6.1.2 Descrizione dell’esempio

In questo esempio si considera una semplice equazione di Poisson.Il dominio e il quadrato (0, 1)× (0, 1).Il tensore e ovunque pari al tensore identita, mentre il termine di sorgentedato da

f = 2π2cos(πx)cos(πy)

Sono date le seguenti condizioni al bordo

p · ν = 0 per (x, y) ∈ 0, 1 × (0, 1)

p = cos(πx) per (x, y) ∈ (0, 1)× 0

p = −cos(πx) per (x, y) ∈ (0, 1)× 1

Si puo mostrare che la soluzione e

p = cos(πx)cos(πy)

6.1.3 Esecuzione e analisi dei risultati

Per compilare e eseguire il file esempio24:

1. Creare una cartella in cui salvare i file necessari alla compilazionedell’esempio.

2. Accedere alla cartella e decomprimere l’archivio ESEMPIO.tar.gz .

3. Se necessario, modificare la variabile Library Path e settarla al per-corso in cui sono stati installate le librerie.

4. Eseguire il comando [$ cmake ./]

5. Eseguire il comando [$ make ]

6. Eseguire il comando [$ ./Esempio]

24Per visualizzarne il codice si veda l’Appendice D.

29

In questo esempio il dominio di calcolo e ottenuto suddividendo prima il do-minio in n step2 quadratini25 di lato h = 1

n step. Ogni quadratino e quindi

suddiviso in 2 triangoli unendo due vertici opposti tramite la diagonale.Per misurare la precisione dell’algoritmo e stato implementata nel file d’e-sempio una versione discreta della norma 2: si valuta

√

∑

(p(xi)− pi)2 neinodi della mesh, dove p(xi) e la soluzione esatta e pi e la soluzione approssi-mata nei suddetti nodi.Per osservare come migliora l’approssimazione, il test viene effettuato sumesh di calcolo sempre piu fitte.Qui sotto vengono riportati i risultati.n step Errore2 0.2859434 0.1667098 0.074515916 0.033679232 0.015831664 0.00765533128 0.00376008256 0.00186478

Per avere le idee piu chiare su cosa sta succedendo, i risultati ottenuti ven-gono plottati su un grafico in scala logaritmica e comparati con una retta dellaforma y = 2x. Come si puo facilmente notare dalla figura 6.1, il metodo e

−6 −5 −4 −3 −2 −1 0

−6

−5

−4

−3

−2

−1

0

x

x

Figura 6.1: Grafico in scala logaritmica dell’errore. Gli asterischi rappresentano i risultati in scalalogaritmica, mentre la linea continua e il grafico della funzione y = 2x.

convergente del prim’ordine.

25Si faccia riferimento al codice riportato in Appendice D.

30

6.1.4 Visualizzazione grafica dei risultati

Come gia accennato, per visualizzare i risultati26 si fara uso del softwareParaview.

1. Avviare paraview [$ paraview].

2. Aprire il file .vtk di output e premere Apply.

Figura 6.2: Campo delle pressioni

3. Per visualizzare “meglio” il campo di pressioni, applicare il filtro Warp

by Scalar.

(a) Soluzione Warped

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1-1

-0.5

0

0.5

1

(b) Soluzione esatta plottatacon Octave

4. Per visualizzare i vettori velocita di Darcy, applicare inizialmente ilfiltro Cell Centers, quindi il filtro Glyph.Si osservi che, come ci si aspetta, il campo di velocita parte da pressionimaggiori per arrivare a pressioni minori (Fig. 6.3).

6.2 Altri esempi

In questa sezione verranno presentati altri due brevi esempi di utilizzo chesi possono trovare rispettivamente in ESEMPIO2.tar.gz e ESEMPIO3.tar.gz.

26L’esempio produrra diversi file di output corrispondenti ai diversi gradi di precisione.Qui si mostrera il procedimento con il file “SoluzioneCHEN 8.vtk”.

31

Figura 6.3: Velocita di Darcy

Le istruzioni per compilarli sono le medesime dell’esempio della sezione prece-dente.

6.2.1 Esempio a dominio generico

Qui di seguito viene presentato un esempio27 in cui il dominio non e unrettangolo (o quadrato). Il file di input e stato generato tramite uno scriptMatlab che si puo trovare in ESEMPIO2.tar.gz.

210.8 210.9 211 211.1 211.2 211.3 211.4 211.5 211.6 211.7 211.8−48.45

−48.4

−48.35

−48.3

−48.25

−48.2

−48.15

−48.1

−48.05

−48

−47.95

Figura 6.4: Mesh di calcolo

In figura (6.4) e riportato il dominio in questione.

27Tratto da [Lie].

32

Il problema che si vuole risolvere e il seguente:

−∆p = 0 ∀x ∈ Ω

p(x, y) = y + xy ∀x = (x, y) ∈ ∂Ω

La soluzione esatta e data da

p(x, y) = y + xy ∀x ∈ Ω

Nonostante i valori non siano ben visibili, si puo notare dalle immagini

(a) Soluzione ottenuta

210.9 211 211.1 211.2 211.3 211.4 211.5 211.6

−48.4

−48.35

−48.3

−48.25

−48.2

−48.15

−48.1

−48.05

−48

−1.028

−1.027

−1.026

−1.025

−1.024

−1.023

−1.022

−1.021

−1.02

−1.019

−1.018

x 104

(b) Soluzione esatta

che la distribuzione di pressioni della soluzione numerica e simile a quelladella soluzione esatta. Se si va a esaminare l’errore nel medesimo modo28

dell’esempio precedente, si ottiene un valore pari a err = 0.00466616, un er-rore accettabile considerando che alcune zone del dominio sono approssimatetramite una griglia relativamente poco fitta.

28Cioe con la norma 2 in versione discreta.

33

6.2.2 Esempio con permeabilita non costante

Adesso si analizza una caso in cui la permeabilita del mezzo non e costantelungo il dominio. Per non complicare troppo lo studio si considera un coeffi-ciente di permeabilita scalare (o equivalentemente un tensore diagonale a ele-menti diagonali uguali). Anche in questo caso la mesh di calcolo (Fig.) e gen-erata attraverso uno script Matlab che si puo trovare in ESEMPIO3.tar.gz.Per la generazione della griglia si e fatto uso di una function Matlab chia-mata distmesh29 in grado di, dato un dominio, fornirne una triangolazionedi Delaunay.Il problema e il seguente:

Figura 6.5: Mesh di calcolo

−∇ · (K(x)∇p(x)) = (x+ y)cos(x+ y)− 2xysin(x+ y) ∀x ∈ Ω

p(x, y) = sin(x+ y) ∀x = (x, y) ∈ ∂Ω

con K(x, y) = xy.La soluzione esatta e data da

p(x, y) = sin(x+ y) ∀x = (x, y) ∈ Ω

Come si puo vedere (Fig. 6.6), anche in questo caso i risultati sonosoddisfacenti e vicini alla soluzione esatta. Un’analisi piu approfondita innorma 2 discreta riporta un errore pari a 0.0026773.

29Implementata da P. Persson e G. Strang. Per ulteriori dettagli si faccia riferimentoalla documentazione ufficiale http://persson.berkeley.edu/distmesh

34

(a) Soluzione ottenuta

1.52

2.53

3.54

4.5522.533.544.55

−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

(b) Soluzione esatta

Figura 6.6: Soluzione dell’esempio a permeabilita variabile nello spazio.

35

Capitolo 7

Conclusione

Come si puo evincere dagli esempi studiati nei precedenti paragrafi, l’algo-ritmo creato e in grado di risolvere con un buon grado di approssimazionegran parte dei problemi della tipologia presentata nel Capitolo 2.Purtroppo, non sempre le cose vanno bene. Un esempio e dato dal seguenteproblema

−∆p = δ(x− x0) ∀x ∈ Ω

p = 0 ∀x ∈ ∂Ω

con Ω = x : |x| ≤ 1.Il modo piu naturale di approssimare il dominio e tramite una griglia circo-lare (Fig. 7.1).

Figura 7.1: Griglia circolare

Nonostante la regolarita della mesh, il metodo CVFE ha particolari difficolta

36

a gestire griglie di questo tipo. Questo perche il CVFE e un metodo chesi basa sulla discretizzazione tramite funzioni polinomiali di primo grado.Per questa ragione, e necessario utilizzare algoritmi costruiti su funzioni piuflessibili. Un esempio e gia stato dato nel Capitolo 3. Infatti il metodo CVFAcorrisponde a questa caratteristica e risulta quindi piu adatto alla risoluzionedi questo particolare problema.Da questa semplice considerazione, si capisce che questo lavoro si presta amolte estensioni. Se ne citeranno solo alcune.La piu naturale evoluzione del codice e l’estensione al caso tridimensionale.Sono gia stati presentati alcuni accorgimenti che lo permettono. Un’idea perampliare il codice puo essere quella di creare una classe NodoCVFE3D cheeredita dalla classe NodoCVFE e aggiunge ad essa un campo che tiene contodella terza componente.Il secondo modo per estendere il programma e, per esempio, considerare icasi non stazionari. Per fare cio, occorrera reimplementare il costruttore e imetodi load() e solve() della classe MeshCVFE.

37

Appendice A

Analisi dei coefficienti ditrasmissibilita

A.1 Principio del massimo e condizioni suffi-

cienti

Per avere le idee piu chiare sulle condizioni da imporre ai coefficienti Tij,bisogna fare, per un momento, un passo indietro e considerare il problemanella formulazione originale. Per semplicita si pongano la permeabilita parial tensore identita e condizioni di Dirichlet su tutta la frontiera ∂Ω ∈ R

2, chesi considera di classe C1 o al piu Lipschitziana30.

−∆u = f in Ω

u = g su ∂Ω

Questo problema ha un’importante proprieta:

Teorema A.1.1. (Principio del massimo)31

Sia u(·) una funzione sufficientemente regolare su Ω ∈ R2. Si ha che

−∆u ≤ 0 ⇒ maxx∈Ω

u(x) = maxx∈∂Ω

u(x)

In particolare se u(·) e negativa sulla frontiera allora e negativa su tutto ildominio.

30Una funzione u(·) si dice Lipschitziana con costante L se

|u(x)− u(y)| ≤ |x− y| ∀x,y ∈ Ω

31Per una dimostrazione si puo consultare un qualsiasi libro di Equazioni alle Derivate

Parziali. Nel caso f = 0 si puo vedere [Salsa].

38

Si osservi che il teorema vale con le disuguaglianze ribaltate e sostituendo laparola positiva alla parola negativa.Sotto particolari condizioni che verranno approfondite piu avanti, questo teo-rema vale anche per funzioni discretizzate32 tramite funzioni lineari a tratticome e stato fatto nella costruzione del metodo oggetto di questa trattazione.Come risultato della discretizzazione, il problema differenziale si riduce allarisoluzione di un sistema lineare

Au = b

dove il vettore u rappresenta la soluzione approssimata nei nodi della grigliadi calcolo.Analizzando la formula (4.9), si possono ottenere le seguenti uguaglianze

Aii =∑

j∈Ωi

Tij

Aij = −Tij

Se A e non-singolare33 e, inoltre, soddisfa

Aii > 0

n∑

j=1

Aij ≥ 0

Aij ≤ 0

Si osserva che queste tre condizioni sono soddisfatte se Tij ≥ 0.Si puo dimostrare che se cio vale, allora A e una M-matrice, cioe invertibilee la sua inversa ha elementi tutti non-negativi.Ma, risultato piu importante, si puo dimostrare che se A e M-matrice allorail problema differenziale discreto associato ad essa rispetta il Principio dimassimo discreto. Tutte queste considerazioni e relative dimostrazioni sipossono trovare in [Ciarlet].Da tutto cio si deduce che il vincolo Tij ≥ 0 e condizione sufficiente affinchesia rispettato il Principio di massimo. Si sottolinea che la condizione e solosufficiente, infatti esistono casi34 in cui condizioni piu deboli garantiscono ilprincipio35.

32Il risultato e noto come Principio del massimo discreto.33Invertibile o , equivalentemente, a determinante non nullo.34Si veda [Korotov].35Esse sono comunque delle condizioni sufficienti. Ne consegue che non e detto che si

ottenga una soluzione sbagliata se non le si rispettano. Ovviamente si cerca sempre diseguirle in quanto nelle applicazioni non si puo lasciare nulla al caso!

39

A.2 Interpretazione geometrica

Il caso bidimensionale e particolarmente semplice e intuitivo poiche am-mette una interpretazione geometrica delle condizioni esposte nella prece-dente sezione.Anche in questo caso si consideri per semplicita la permabilita pari al tensoreidentita.

Figura A.1: Elementi che contribuiscono a un generico Tij .

Riprendendo le formule (4.5) e (4.8) e considerando un singolo elemento dellamesh, si puo provare che vale l’uguaglianza

Tij = −aiaj + bibj

4|M |

cioe (Fig. A.1)

Tij =|mkmj||mkmi|cos(θk)

4|M |=

cot(θk)

2

Ora, considerando il fatto che ad ogni Tij contribuiscono due elementi Mdifferenti della mesh, si ha

Tij =cot(θ1) + cot(θ2)

2

e la condizione Tij ≥ 0 si traduce nella condizione

θ1 + θ2 < π

Nella condizione di tensore di permeabilita costante, questo vincolo equivalea richiedere che il dominio di calcolo sia una triangolazione di Delaunay,cioe ogni elemento (triangolare) della griglia e tale che il cerchio ad essocircoscritto non contiene altri nodi della mesh.

40

Appendice B

Eccezioni

Qui si riportano tutte le eccezioni che potrebbero essere lanciate dalla libreria.

1 namespace MeshCvfe2 namespace ErroriCVFE3 s t r u c t ErroreASCII : pub l i c except ion 4 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

Formato f i l e non r i c o n o s c i u t o . I l f i l e deve e s s e r e inASCII ! ” ;

5 ;6 s t r u c t ErroreElemAdiacent i : pub l i c except ion 7 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

Troppi e l ement i ad i a c en t i . Cont ro l l a r e i l f i l e ! ” ;8 ;9 s t r u c t Er r o r eF i l e : pub l i c except ion

10 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Erroren e l l ’ apertura de l f i l e ! ” ;

11 ;12 s t r u c t ErroreCel lTypes : pub l i c except ion 13 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

numero d i da t i CELL TYPES non cor r i spondente a l numeroprecendentemente d i c h i a r a t o d i CELLS” ;

14 ;15 s t r u c t ErrorePointData : pub l i c except ion 16 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

numero d i POINT DATA non cor r i spondente a l numeroprecendentemente d i c h i a r a t o d i POINTS” ;

17 ;18 s t r u c t ErroreCel lData : pub l i c except ion 19 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

numeri d i CELL DATA non co r r i s p onden t i a l numeroprecedentemente d i c h i a r a t o d i CELLS” ;

20 ;21 s t r u c t ErroreElemento : pub l i c except ion

22 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :Deve e s s e r e un t r i a n g o l o (CELL TYPE=5)” ;

23 ;24 s t r u c t ErroreLat i : pub l i c except ion 25 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

numero d i l a t i e r r a t o . La c e l l a deve e s s e r e unt r i a n g o l o ! ” ;

26 ;27 s t r u c t ErroreGrid : pub l i c except ion 28 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

datase t non g e s t i b i l e . I da t i devono e s s e r e o r g an i z z a t iin una g r i g l i a non s t r u t t u r a t a (UNSTRUCTUREDGRID) ” ;

29 ;30 s t r u c t ErroreNaN : pub l i c except ion 31 v i r t u a l const char ∗ what ( ) const throw ( ) r e turn ”Errore :

s t a t o generato un Not−A−Number” ;32 ;33 34

42

Appendice C

La classe Point2D

Qui viene presentata la classe Point2D (da cui eredita NodoCVFE). Per faruso della classe bisogna includere l’header Point2D.hpp .La classe e all’internodel namespace POINT2D. Per poter far uso della relativa libreria occorre, almomento della compilazione, aggiungere l’opzione [$ -lPoint2D ]

1 c l a s s Point2D2 pub l i c :3 Point2D ( double X, double Y) : x (X) , y (Y) 4 Point2D ( )=de f au l t ;5 double coordx ( ) const ;6 double coordy ( ) const ;7 void se tx ( double ) ;8 void se ty ( double ) ;9 double d i s tanza ( const Point2D&) const ;

10 f r i e nd ostream& operator<<(ostream&, const Point2D ) ;11 ˜Point2D ( )=de f au l t ;12 pr i va t e :13 double x=0;14 double y=0;15 ;

Data la semplicita della classe non si rendono necessarie ulteriori spiegazioni.

43

Appendice D

Esempio

Qui si riporta il codice dell’esempio.

1 #inc lude <MeshCVFE. hpp>2 #inc lude <Point2D . hpp>3 #inc lude <f stream>

4 #inc lude <sstream>

5 #inc lude <s t r i ng>6 #inc lude <math . h>7 us ing namespace POINT2D;8 us ing namespace MeshCvfe ;9 us ing namespace std ;

10 void GeneraEsempioCHEN( in t n s t ep ) ;11

12 i n t main ( ) 13 i n t num step=1;14 MeshCVFE mesh ;15 ofstream l o g f i l e ;16 //Apertura l o g f i l e17 l o g f i l e . open ( ”ExampleCHEN . log ” , i o s : : out | i o s : : trunc ) ;18 l o g f i l e <<”Log F i l e d e l l ’ esempio . Di s e gu i t o sono r i p o r t a t i i l

numero d i pas s i , l ’ e r r o r e e i l tempo di e s e cuz i one . Buonal e t t u r a ! ”<<endl ;

19 f o r ( i n t e=1;e<9; e++)20 num step∗=2;21 l o g f i l e <<”Numero d i p a s s i : ”<<num step<<”\ t \ t \ t ” ;22 //Creaz ione d e l l ’ esempio con numero d i p a s s i pa r i a 2ˆ e23 GeneraEsempioCHEN( num step ) ;24 //Creaz ione e r i s o l u z i o n e d e l l a mesh25 mesh . load ( ”EsempioCHEN . vtk” ) ;26 mesh . s o l v e ( ) ;27 // Calco lo e s c r i t t u r a d e l l ’ e r r o r e28 double max norm=0;29 f o r ( i n t i =0; i<mesh . showNumNodi ( ) ; i++)30 double t ;

31 double x , y ;32 x=mesh [ i ] . coordx ( ) ;33 y=mesh [ i ] . coordy ( ) ;34 t=fabs (mesh [ i ] . p r e s su r e ( )−cos (M PI∗x ) ∗ cos (M PI∗y ) ) ;35 i f (max norm<t ) max norm=t ;36 37 l o g f i l e <<”Errore : ”<<max norm<<endl ;38 s t r i ng s t r eam output ;39 // Sa lva tagg io d e l l a s o l u z i on e40 output<<”SoluzioneCHEN ”<<e<<” . vtk” ;41 mesh . pr intVtk ( output . s t r ( ) ) ;42 43 r e turn 0 ;44

Il file input per l’esempio e generato grazie a questo algoritmo.

1 #inc lude <iostream>

2 #inc lude <f stream>

3 #inc lude <c s t r i ng>4 #inc lude <ErroriCVFE . hpp>5 #inc lude <math . h>6 #inc lude <Point2D . hpp>7

8

9 void GeneraEsempioCHEN( in t n s t ep ) 10 //Apertura f i l e11 us ing namespace POINT2D;12 us ing namespace MeshCvfe ;13 double h=1/(double ) n s t ep ;14 typede f i n t num type ;15 ofstream ou t p u t f i l e ;16 o u t p u t f i l e . open ( ”EsempioCHEN . vtk” , i o s : : out | i o s : : trunc ) ;17 i f ( o u t p u t f i l e . good ( ) ) 18 // I n t e s t a z i o n e19 ou t pu t f i l e<<”# vtk DataFi le Vers ion 3 .0\n” ;20 // Des c r i z i one21 ou t pu t f i l e<<”Esempio preso dal l i b r o CHEN−CVFE ”<<endl ;22 //Formato23 ou t pu t f i l e<<”ASCII\n\n” ;24 // St rut tura25 ou t pu t f i l e<<”DATASET UNSTRUCTUREDGRID\n” ;26 //Coordinate nodi27 num type num nodi=(n s tep+1)∗( n s t ep+1) ;28 ou t pu t f i l e<<”POINTS ”<<num nodi<<” f l o a t \n” ;29 vector<Point2D> v ps ;30 v ps . r e s e r v e ( num nodi ) ;31 f o r ( num type i =0; i<num nodi ; i++)

45

32 double x , y ;33 x=( i%(n s tep+1) ) ∗h ;34 y=1−(( i n t ) ( i /( n s t ep+1) ) ) ∗h ;35 ou t pu t f i l e<<x<<” ”<<y<<” 0\n” ;36 Point2D p(x , y ) ;37 v ps . push back (p) ;38 39 // Ve r t i c i Ce l l e ( t r i a n g o l i )40 num type num elem=n step ∗ n s tep ∗2 ;41 ou t pu t f i l e<<”\nCELLS ”<<num elem<<” ”<<num elem∗4<<endl ;42 vector<double> f i n c e n t r i ;43 f i n c e n t r i . r e s e r v e ( num elem ) ;44 f o r ( num type i =0; i<n s tep ; i++)45 f o r ( num type j =0; j<n s tep ; j++)46 i n t v1 , v2 , v3 ;47 v1=j+(n s tep+1)∗ i ;48 v2=j+(n s tep+1)∗( i +1)+1;49 v3=j+(n s tep+1)∗ i +1;50 double a , b , c , p ;51 a=sq r t ( ( v ps [ v1 ] . coordx ( )−v ps [ v2 ] . coordx ( ) ) ∗( v ps [ v1 ] .

coordx ( )−v ps [ v2 ] . coordx ( ) )+(v ps [ v1 ] . coordy ( )−v ps [v2 ] . coordy ( ) ) ∗( v ps [ v1 ] . coordy ( )−v ps [ v2 ] . coordy ( ) ) ) ;

52 b=fabs ( v ps [ v3 ] . coordy ( )−v ps [ v2 ] . coordy ( ) ) ;53 c=fabs ( v ps [ v3 ] . coordx ( )−v ps [ v1 ] . coordx ( ) ) ;54 p=(a+b+c ) ;55 double x , y ;56 x=(a∗ v ps [ v3 ] . coordx ( )+b∗ v ps [ v1 ] . coordx ( )+c∗ v ps [ v2 ] .

coordx ( ) ) /p ;57 y=(a∗ v ps [ v3 ] . coordy ( )+b∗ v ps [ v1 ] . coordy ( )+c∗ v ps [ v2 ] .

coordy ( ) ) /p ;58 f i n c e n t r i . push back (2∗ (M PI∗M PI) ∗ cos (M PI∗x ) ∗ cos (M PI∗

y ) ) ;59 ou t pu t f i l e<<”3 ”<<v1<<” ”<<v2<<” ”<<v3<<endl ;60 v1=j+(n s tep+1)∗ i ;61 v2=j+(n s tep+1)∗( i +1) ;62 v3=j+(n s tep+1)∗( i +1)+1;63 a=sq r t ( ( v ps [ v1 ] . coordx ( )−v ps [ v3 ] . coordx ( ) ) ∗( v ps [ v1 ] .

coordx ( )−v ps [ v3 ] . coordx ( ) )+(v ps [ v1 ] . coordy ( )−v ps [v3 ] . coordy ( ) ) ∗( v ps [ v1 ] . coordy ( )−v ps [ v3 ] . coordy ( ) ) ) ;

64 b=fabs ( v ps [ v1 ] . coordy ( )−v ps [ v2 ] . coordy ( ) ) ;65 c=fabs ( v ps [ v2 ] . coordx ( )−v ps [ v3 ] . coordx ( ) ) ;66 p=(a+b+c ) ;67 x=(a∗ v ps [ v2 ] . coordx ( )+b∗ v ps [ v3 ] . coordx ( )+c∗ v ps [ v1 ] .

coordx ( ) ) /p ;68 y=(a∗ v ps [ v2 ] . coordy ( )+b∗ v ps [ v3 ] . coordy ( )+c∗ v ps [ v1 ] .

coordy ( ) ) /p ;69 f i n c e n t r i . push back (2∗ (M PI∗M PI) ∗ cos (M PI∗x ) ∗ cos (M PI∗

y ) ) ;70 ou t pu t f i l e<<”3 ”<<v1<<” ”<<v2<<” ”<<v3<<endl ;

46

71 72 73 ou t pu t f i l e<<endl ;74 // Tipo log ia c e l l e75 ou t pu t f i l e<<”CELL TYPES ”<<num elem<<endl ;76 f o r ( num type i =0; i<num elem ; i++)77 ou t pu t f i l e<<”5\n” ;78 79

80 // P r e s s i o n i ne i nodi81 ou t pu t f i l e<<endl<<”POINT DATA ”<<num nodi<<endl ;82 ou t pu t f i l e<<”SCALARS incogn i t a i n t \n” ;83 ou t pu t f i l e<<”LOOKUPTABLE de f au l t \n” ;84 f o r ( num type i =0; i<n s tep+1; i++)85 f o r ( num type j =0; j<n s tep+1; j++)86 i f ( i==0 | | i==n step ) 87 ou t pu t f i l e<<”0”<<endl ;88 e l s e 89 ou t pu t f i l e<<”1”<<endl ;90 91 92 93 ou t pu t f i l e<<endl ;94 ou t pu t f i l e<<”SCALARS pr e s s i on e f l o a t \n” ;95 ou t pu t f i l e<<”LOOKUPTABLE de f au l t \n” ;96 f o r ( num type i =0; i<n s tep+1; i++)97 f o r ( num type j =0; j<n s tep+1; j++)98 i f ( i==0)99 ou t pu t f i l e<<−cos (M PI∗ j ∗h)<<endl ;

100 e l s e 101 i f ( i==n step ) 102 ou t pu t f i l e<<cos (M PI∗ j ∗h)<<endl ;103 e l s e 104 ou t pu t f i l e <<0<<endl ;105 106 107 108 109 // Sorgente110 ou t pu t f i l e<<endl<<”CELL DATA ”<<num elem<<endl ;111 ou t pu t f i l e<<”SCALARS sorgente f l o a t \n” ;112 ou t pu t f i l e<<”LOOKUPTABLE de f au l t \n” ;113 f o r ( num type i =0; i<num elem ; i++)114 ou t pu t f i l e<< f i n c e n t r i [ i ]<<endl ;115 116 //Tensor i117 ou t pu t f i l e<<”TENSORS permeab i l i t y f l o a t \n” ;118 f o r ( num type i =0; i<num elem ; i++)119 ou t pu t f i l e<<”1 0 0”<<endl ;

47

120 ou t pu t f i l e<<”0 1 0”<<endl ;121 ou t pu t f i l e<<”0 0 1”<<endl<<endl ;122 123 o u t p u t f i l e . c l o s e ( ) ;124 r e turn ;125 e l s e throw ErroriCVFE : : E r r o r eF i l e ( ) ;126

48

Bibliografia

[Chavent] G. Chavent, J. Jaffre, Mathematical Models and Finite Elements

For Reservoir Simulation, North Holland, 1986

[Chen] Z. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase

Flows in Porous Media, Siam, 2006

[Ciarlet] P.G. Ciarlet, Discrete maximum principle for finite difference

operators (in Aequationes Math. 4 (1970), 338-352)

[Fanchi] J. Fanchi, Principles of Applied Reservoir Simulation, 2 ed., GulfProfessional Publishing, 2001

[Lie] K.A.Lie, S.Krogstad, I.S. Ligaarden, J.R. Natvig, H.M. Nilsen,B. Skaflestad, Open Source MATLAB Implementation of Consistent

Discretisations on Complex Grids

[Lippman] S.B. Lippman,J. Lajoie, B.E.Moo, C++ Primer, 5th ed.,Addison-Wesley,2012

[Korotov] S. Korotov, M. Krizek, P. Neittaanmaki, Weakened acute type

condition for tetrahedral triangulations and the discrete maximum

principle

[Pinder] G. Pinder, W. Gray, Essentials of Multiphase Flow and Transport

in Porous Media, wiley,2008

[Quarteroni] A. Quarteroni, Modellistica numerica per problemi differenziali,4 ed., Springer, 2008

[Salsa] S. Salsa, Equazioni a derivate parziali, 2 ed., Springer, 2010

49