Volúmenes Finitos en varias dimensiones

Métodos Numéricos Avanzados. 2023-2024

Author
Affiliation

David Gómez-Castro

Universidad Complutense de Madrid

Problemas en dimensión superior

Diviendo el espacio en regiones

Hasta ahora hemos dividido un intervalor \(\mathbb R\) mediante una partición.

Si vamos a trabajar con un problema un un sub-dominio de \(\mathbb R^d\) el primer paso es partirlo en regiones pequeñas, debemos extender esta idea a “mallas”.

Una malla es una división de nuestro dominios en polígonos convexos.

Dado un dominio en dimensión \(d \ge 2\) construir una descomposición es complicado. Hay diferentes software que resuelven este problema

Ejemplo de Gmsh

Supuesto que ya tenemos una malla, la pregunta es ¿cómo escribo mis ecuaciones sobre la malla?

Malla

Seguimos la presentación de (Eymard, Gallouët & Herbin, 1997).

Una malla admisible de un conjunto \(\Omega\) es un par \((\mathcal T, \mathcal E)\) donde \(\mathcal T\) es es una familia de polígonos conexos \(K\) (llamados volúmenes de control), \(\mathcal E\) una familia de subconjuntos \(\sigma\) de \(\overline \Omega\) contenidos en hiperplanos de \(\mathbb R^d\) (llamados aristas o lados) tales que:

  • Son disjuntos: \(K , L \in \mathcal T\) distintos, entonces \(K \cap L = \emptyset\).

  • Para cada \(K\) existe \(\mathcal E_K \subset \mathcal E\) con \(\partial K = \cup_{\sigma \in \mathcal E_K} \overline \sigma\).

  • La unión de sus cierres genera \(\mathcal T\): \(\bigcup_{K \in \mathcal T} \overline K = \overline \Omega\)

  • Entre cualesquiera \(K, L \in \mathcal L\) distintos la “interfaz” \(\overline K \cap \overline L\) tiene medida \(\mathcal H^{d-1}\) igual 0 o es \(\overline \sigma\) con \(\sigma \in \mathcal E\). Llamamos a \(\sigma = K|L\).

Llamamos \(\mathcal N(K)\) el conjunto de \(L\) tales que \(\mathcal H^{d-1} (K|L) \ne 0\).

Además, \(\mathcal E_{\textrm{ext}} = \{ \sigma \in \mathcal E : \sigma \subset \partial \Omega \}\).

Llamamos \(\mathbf n_{K,\sigma}\) el vector normal a \(\sigma\), exterior a \(K\).

Problema continuo

Sobre un volumen de control

\[ \frac{\diff}{\diff t} \int_{K} \rho(t,x) \diff x = - \sum_{\sigma \in \mathcal E_K} \int_{\sigma} \bm \flux(t,x) \cdot \bm n_{K,L}(x) \diff S(x) . \]

Esto significa que puedo buscar ecuaciones para \[ \rho_K (t) \approx \frac{1}{|K|} \int_{K} \rho(t,x) \diff x. \]

Y para cada polígono adyacente \(L\) hemos de construir \[ F_{K,\sigma} \approx \bm \flux(t,x) \cdot \bm n_{K,\sigma} (x) \]

La conservación de masa entre celdas corresponde a \(F_{L, K|L} = - F_{K, K|L}\).

La condición de flujo nulo es \(F_{K,\sigma} = 0\) if \(\sigma \in \mathcal E_{\textrm{ext}}\).

Podemos plantear un problema discreto \[ m(K) \frac{\diff}{\diff t} \rho_K = - \sum_{\sigma \in \mathcal E_K} m(\sigma) F_{K,\sigma} \]

donde \(m(K) = \mathcal H^d (K)\) y \(m(K|L) = \mathcal H^{d-1} (K|L)\)

Mallas rectangulares

Si trabajamos en un dominio cuadrado, \([-1,1]^2\) por ejemplo, podemos construir los polígonos \[ K_{ij} = [ x_{i-\frac 1 2} , x_{i+\frac 1 2} ] \times [y_{j-\frac 1 2} , y_{j+\frac 1 2}]. \]

Escribiendo \(\bm \flux = (G,H)\) entonces \[ \frac{\diff \rho_{ij}}{\diff t} (t) = - \frac{G_{i+\frac 1 2 , j} (t) - G_{i-\frac 1 2, j} (t)}{\Delta x_i} - \frac{H_{i,j+\frac 1 2}(t) - H_{i,j-\frac 1 2} (t)}{\Delta y_j} \]

Seguimos pudiendo aplicar argumentos de up-winding y monotonía.

Experimentos numéricos

Tomamos \(y = x\).

Code 
using Plots, Printf
function create_mesh(N)
    x  = LinRange(-1.0,1.0,N+1); x = [x[1:end-1];]

    ρ0 = zeros(N,N)
    for i in 1:N, j in 1:N
        ρ0[i,j] = max(0.0, 1.1*exp(-(x[i]^2 + x[j]^2)/.05) - .1 )
    end

    Δx = x[2] - x[1] 

    return (x,Δx,ρ0)
end

function myplot(x,ρ,t,thetitle)
    surface(x,x,ρ,
        title= thetitle*"\n"*@sprintf("t = %0.3f, max = %0.3f, mass = %0.3f", 
            t,maximum(ρ),(x[2]-x[1])^2*sum(ρ)),
        colorbar=false,zlim=[0,1],
        xlabel="x",ylabel="y") 
end;
N = 32; (x,Δx,ρ) = create_mesh(N)

myplot(x,ρ,0,"Dato inicial")

Transporte a velocidad constante y Euler explícito

Pongamos que \(\bm \flux = (a_1, a_2) \rho\).

Podemos hacer up-winding en cada una de las componentes: \[ \begin{aligned} G_{i+\frac 1 2, j} &= (a_1)_+ \rho_{ij} + (a_1)_- \rho_{i+1,j}, \\ H_{i, j+\frac 1 2} &= (a_2)_+ \rho_{ij} + (a_2)_- \rho_{i,j+1} \end{aligned} \]

Método de Euler explícito con up-winding

Code 
N = 32; (x,Δx,ρ) = create_mesh(N)
a1  = 1; a2  = 0.25;
Δt = Δx / (a1 + a2) # CFL

periodic(i) = mod(i-1,N)+1

T = 2
@gif for n=1:ceil(T/Δt) 
    ρn = copy(ρ)
    myplot(x,ρ,(n-1)*Δt,"Transporte. Euler explícito")

    for i in 1:N, j in 1:N
        ρ[i,j] = ρn[i,j] - 
            a1 * Δt/Δx * (ρn[i,j] - ρn[periodic(i-1),j] ) -
            a2 * Δt/Δx * (ρn[i,j] - ρn[i,periodic(j-1)] )
    end

end

Ecuación de calor

El flujo sigue siendo \[ \bm \flux = -\nabla \rho = -\left( \frac{\partial \rho}{\partial x} , \frac{\partial \rho}{\partial y} \right) \]

Como en la malla cuadrada las normales son \((1,0)\) y \((0,1)\) (salvo el signo), es fácil hacer los cálculos. Hacemos las diferencias centradas \[ \frac{\diff}{\diff t} \rho_{ij} (t) = - \frac{2 \rho_{ij} (t) - \rho_{i+1,j}(t) - \rho_{i-1,j}(t) }{(\Delta x)^2} - \frac{2 \rho_{ij} (t) - \rho_{i,j+1}(t) - \rho_{i,j-1}(t) }{(\Delta y)^2} \]

El método de Euler explícito nos lleva a una condición \[ 2\frac{\Delta t}{(\Delta x)^2} + 2\frac{\Delta t}{(\Delta y)^2} \le 1. \] Dado que ahora trabajamos con matrices \(N \times N\) es computacionalmente muy costosa.

Ecuación de Calor. Método de Euler implícito

Querríamos usar el método implícito

\[ \frac{\rho_{ij}^{n+1} - \rho_{ij}^n}{\Delta t} = - \frac{2 \rho_{ij}^{n+1} - \rho_{i+1,j}^{n+1} - \rho_{i-1,j}^{n+1} }{(\Delta x)^2} - \frac{2 \rho_{ij}^{n+1} - \rho_{i,j+1}^{n+1} - \rho_{i,j-1}^{n+1} }{(\Delta y)^2} \]

En dimensión 1, era fácil despejar \(A \bm \rho^{n+1} = \bm \rho^n\).

Pero ahora \(\rho^n\) es una matriz, y ya no es tan fácil hacer esta escritura.

Producto de Kronecker

La manera natural de pensar en \(\rho_{ij}\) es como matriz.

Sin embargo, las aplicaciones sobre matrices son tensores de orden 3. Son difíciles de manejar.

Para evitar esto, hacemos un re-shape de la matriz \(N \times M\) a vector \(NM\) \[ \begin{aligned} \mathrm{vec} : \mathcal M_{N,M} &\longrightarrow \mathbb R^{NM} \\ (\rho_{ij}) &\longmapsto \Big( (\rho_{11}, \rho_{12}, \cdots, \rho_{1M}), (\rho_{21}, \rho_{22}, \cdots, \rho_{2M}), \cdots , (\rho_{N1}, \rho_{N2}, \cdots, \rho_{NM}) \Big)^t \end{aligned} \]

Podemos pensar en esto es tomar ver esta operación como el “desmontaje” de un producto tensorial. Sean \(u_i\) la base canónica de \(\mathbb R^N\) y \(v_j\) la de \(\mathbb R^M\). Asociando el producto tensorial \(u_i \times v_j\) con la matriz que tiene sólo un 1 y en la entrada (\(i,j\)), podemos escribir \[ \rho = \sum_{i,j} \rho_{ij} u_i \otimes v_j. \]

Así \[ \mathrm{vec} ( u_i \otimes v_j ) = e_{M(i-1) + j}. \]

Producto tensorial y matrices por bloques

Dado que tenemos un operador (“diferencial”) operando sólo en \(u_i\) querríamos ver cómo es la matriz que corresponde a aplicarlos los dos. Querremos hacer lo mismo operando sólo en \(v_j\). Dados operadores \(A\) y \(B\), se define el operador \(A \otimes B\) de la siguiente forma \[ (A \otimes B) \rho = \sum_{i,j} \rho_{ij} (A u_i) \otimes (B v_j). \] Si queremos expresar esto como una matriz \(C \in M_{NM, NM}\), la escribimos por bloques \[ C = \begin{pmatrix} C_{11} & \cdots & C_{1N} \\ \vdots & \ddots & \vdots \\ C_{N1} & \cdots & C_{NN} \end{pmatrix}, \qquad C_{ij } \in \mathcal M_{M,M} \] Por un lado \[ (A\otimes B) (u_k \otimes v_\ell) = \left( \sum_i a_{ik} u_i \right) \otimes \left( \sum_j b_{j\ell} v_j \right) = \sum_{i,j} a_{ik}b_{j \ell} (u_i \otimes v_j) . \] Comprobamos la entrada \(K = M(k-1) + \ell\)-ésima fila \[ (C \mathrm{vec} (u_k \otimes v_\ell) )_K = ( C e_{M(k-1) + \ell} )_K = C_{K , M(k-1) + \ell } = (C_{ik})_{j \ell}, \] de modo que \[ C_{ik} = a_{ik} B. \] A la operación \(\otimes : \mathcal M_{N} \times \mathcal M_M \to \mathcal M_{NM}\) \[ A \otimes B = \begin{pmatrix} a_{11}B & \cdots & a_{1N} B \\ \vdots & \ddots & \vdots \\ a_{N1} B & \cdots & a_{NN} B \end{pmatrix} \] se la denomina producto de Kronecker de las matrices.

Observación. Ocurre que si \(X\) es una matriz entonces

\[ (A \otimes B) \mathrm{vec} (X) = \mathrm{vec}( B^t X A ) \]

Programando \(\mathrm{vec}\) y \(\mathrm{vec}^{-1}\)

function mat2vec(A)
    (N,M) = size(A)
    u = zeros(N*M,1)
    for i in 1:N, j in 1:M
        u[M*(i-1)+j] = A[i,j]
    end
    return u
end

function vec2mat(u,N)
    M = div(length(u),N)
    A = zeros(N,M)
    for i in 1:N, j in 1:M
        A[i,j] = u[M*(i-1) + j] 
    end
    return A
end;

Ecuación de transporte. Euler implícito.

Construímos la matriz \(A\) tal que \[A \mathrm{vec}(\rho^{n+1}) = \mathrm{vec}(\rho^n)\] mediante el producto de Kronecker.

Como es muy grande (\(N^2 \times N^2 = N^4\)) pero está llena de 0s, la almacenamos como matriz sparse

N=32; C  = 0.5;

D = zeros(N,N)
for i=1:N   D[i,i]   =  1 end
for i=1:N-1 D[i+1,i] = -1 end
D[1,N] = -1

using LinearAlgebra, Kronecker, SparseArrays
Id  = Matrix(1.0I, N, N)
A   = sparse( Matrix(1.0I, N*N, N*N) 
        + C * ( a1 * D ⊗ Id + a2 * Id ⊗ D ) )
1024×1024 SparseMatrixCSC{Float64, Int64} with 3072 stored entries:
⎡⣷⣅⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠦⎤
⎢⠈⠻⣷⣔⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠈⠻⣷⣄⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠈⠻⣷⣌⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠻⣷⣄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣅⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣔⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣄⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣌⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣔⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣷⣅⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣌⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣅⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣔⠀⡀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣄⠄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣌⢀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣄⠂⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣅⠠⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⣷⣔⎦
ρ = mat2vec(ρ)
T = 3
@gif for n=1:ceil(T/Δt)
    myplot(x,vec2mat(ρ,N),(n-1)*Δt,
        "Transporte. Euler implícito")

    global ρ = A \ ρ
end

Dimensional Splitting

La maldición de la dimensionalidad

Cuando se trabaja en dimensión espacial \(d>1\), especialmente con métodos implícitos, la complejidad computacional aumenta. Tomando \(N\) puntos en cada dimensión, el tamaño de \(\bm \rho\) es \(N^d\). Si se trabaja con un método implícito, las matrices que aparecen son de tamaño \(N^d\times N^d\) (generalmente sparse).

Una forma de lidiar con este problema es utilizar dimensional splitting.

Idea

Empecemos por la idea del splitting. Por establecer un marco de trabajo, consideremos el problema \[ \frac{\partial u}{\partial t} = Au + Bu, \qquad u(0) = u_0, \] donde \(A,B\) son “operadores” lineales o no-lineales en \(u\).

Una idea que surge es resolver iterativamente problemas de la forma \[ \begin{aligned} \frac{\partial v}{\partial t}&= Av, \qquad v(0) = v_0,\\ \frac{\partial w}{\partial t}& = Bw, \qquad w(0) = w_0, \end{aligned} \]

Para ser precisos queremos resolver \[ \begin{aligned} \frac{\partial v}{\partial t}&= Av, \text{ para } t \in [0,1/n] \qquad v(0) = u_0,\\ \frac{\partial w}{\partial t}& = Bw, \text{ para } t \in [0,1/n]\qquad w(0) = v(1/n),\\ \\ \frac{\partial v}{\partial t}&= Av, \text{ para } t \in [1/n,2/n] \qquad v(1/n) = w(1/n),\\ \frac{\partial w}{\partial t}& = Bw, \text{ para } t \in [1/n,2/n] \qquad w(1/n) = v(2/n), \\ \vdots \end{aligned} \]

Con la notación de semigrupos, denotemos \[ u(t) = S(t)u_0, \qquad v(t) = S_A (t) v_0, \qquad \text{y} w(t) = S_B (t) w_0 \] Definamos \[ u_k (t) = S_B (\tfrac t k ) \circ S_A (\tfrac t k ) \circ \overset k \cdots \circ S_B (\tfrac t k ) \circ S_A (\tfrac t k ). \] En ocasiones se cumple que \[ u_k (t) \to u(t) , \qquad \text{ as } k \to \infty. \]

Esta idea se puede extender a muchos operadores. La formulación más elegante del resultado es \[\begin{equation} u(t) = \lim_{k\to \infty} \Big[ S_B (\tfrac t k ) S_A (\tfrac t k ) \Big]^k u_0 \end{equation}\] A esta fórmula se la conoce como fórmula de Trotter-Kato.

La intuición viene del caso lineal, en el que \[ e^{(A+B)t} = \lim_k \left( e^{B \tfrac t k } e^{A\tfrac t k } \right) ^k \] Esta fórmula es exacta para, por ejemplo, matrices.

Dimensional splitting

Si queremos resolver un problema de forma \[ \frac{\partial \rho}{\partial t} = - \frac{\partial }{\partial x} f(\rho) - \frac{\partial }{\partial y} g(\rho) \]

Podemos tomar \[ A u = - \frac{\partial}{\partial x} f(u) , \qquad B u = - \frac{\partial }{\partial y} g(u). \]

y construir el método en dos pasos \[ \begin{aligned} \rho^*_{ij} &= \rho_{ij}^n - \frac{\Delta t}{\Delta x} \left({ f_{i+\frac 1 2, j}^{n+1} - f_{i-\frac 1 2,j}^{n+1} }\right) \\ \rho^{n+1}_{ij} &= \rho_{ij}^* - \frac{\Delta t}{\Delta y} \left( g_{i, j + \frac 1 2}^{n+1} - g_{i , j-\frac 1 2}^{n+1} \right) \end{aligned} \]

Para calcular \(g_{ij}^{n+1}\) utilizamos \(\rho^*\) o \(\rho^{n+1}\), pero no \(\rho^n\).

Si el método es implícito, en casa paso tenemos que resolver:

  • \(M\) sistemas de \(N\) ecuaciones para calcular \(\bm \rho^*\)
  • \(N\) sistemas de \(M\) ecuaciones para calcular \(\bm \rho^{n+1}\)

Si \(N = M\) en el método de Newton hay \(2N\) matrices \(N \times N\). Pero no una matriz \((N^2) \times (N^2)\).

Observación. Strang demostró que la convergencia de \[ u(t) = \lim_{k\to \infty} \Big[ S_A (\tfrac t {2k} ) S_B (\tfrac t k ) S_A (\tfrac t {2k} ) \Big]^k u_0 \] es mejor que la de la fórmula original de Trotter-Kato.

Ecuación de transporte con condiciones periódicas

Construímos las matrices en ambas direcciones

T = 3.0; N=32; (x,Δx,ρ) = create_mesh(N); C  = 0.5; Δt = C*Δx

D = zeros(N,N)
for i=1:N   D[i,i]   =  1 end
for i=1:N-1 D[i+1,i] = -1 end
D[1,N] = -1

Id = Matrix(1.0I, N, N)
A1 = sparse( Id + C * a1 * D )
A2 = sparse( Id + C * a2 * D )

@gif for n=1:ceil(T/Δt)
    global ρ
    myplot(x,ρ, (n-1)*Δt, "Transporte. Euler implícito. Dimensional splitting")
    # Difusión en x con y fija
    for j in 1:N
        ρ[:,j] = A1 \ ρ[:,j] 
    end    
    # Difusión en y con x fija
    for i in 1:N
        ρ[i,:] = A2 \ ρ[i,:] 
    end
end

Ecuación de calor con condiciones periódicas

T = 3e-2; N=32; (x,Δx,ρ) = create_mesh(N); Δt = 1e-3; C  = Δt / (Δx)^2

A = sparse(SymTridiagonal((1+2*C)*ones(N),-C*ones(N-1)))
A[1,N] = A[N,1] = -C;

anim = @animate for n=1:ceil(T/Δt) 
    global ρ
    myplot(x,ρ,(n-1)*Δt,"Calor. Euler implícito. Dimensional splitting")

    # Difusión en x con y fija
    for j in 1:N
        ρ[:,j] = A \ ρ[:,j] 
    end    

    # Difusión en y con x fija
    for i in 1:N
        ρ[i,:] = A \ ρ[i,:] 
    end
end
gif(anim,fps=3) # Reduce la velocidad del gif

Mallas generales

Habitualmente el dominio no es un cuadrado y, en general, no es deseable usar mallas cuadradas.

Sin embargo, en el caso de mallas generales más difícil hacer una teoría uniforme.

Presentamos un ejemplo.

La ecuación de transporte

Adaptado de (Eymard, Gallouët & Herbin, 1997, Sección 25)

Pensemos en la ecuación de transporte, donde \(\mathbf F(t,x) = \rho(x) \mathbf v(x)\) donde \(v\) es conocido.

Si uno quiere hacer una aproximación del flujo a través del borde, por ejemplo

Si \(\sigma = K|L\) la velocidad con up-wind \[ v_{K,\sigma} = \frac{1}{m(\sigma)}\int_{\sigma} (\mathbf v (x) \cdot \mathbf n_{K,\sigma})_+ dS(x) \]

el up-winding en el caso continuo corresponde a

\[ F^n_{K,\sigma} = \rho_K^n v_{K,\sigma} - \rho_{L}^n v_{L,\sigma} \]

La condición de flujo nulo en el exterior es \(F_{K,\sigma}= 0\) para \(\sigma \in \mathcal E_{\textrm{ext}}\)

\[ m(K) \frac{\rho_K^{n+1} - \rho_K^n}{\Delta t} = - \sum_{\sigma \in \mathcal E_K \cap \mathcal E_{\textrm{int}}} m(\sigma) F^n_{K,\sigma} . \]

La ecuación de calor

Adaptado de (Eymard, Gallouët & Herbin, 1997, Capítulo 3)

Para esto debe asumir que nuestra malla es tal que existen \(x_K \in \overline K\) tales que \([x_K,x_L] \perp K|L\) y \(x_K \ne x_L\) si \(K \ne L\).

Así \(\mathbf F = -\nabla \rho\) se puede aproximar \(F_{K,\sigma} \approx -\nabla \rho \cdot n_{K,\sigma}\): \[ F_{K,\sigma} = - \frac{\rho_{L} - \rho_K}{ |x_L - x_K| }, \qquad \sigma = K|L . \]

Añadiendo la condición de flujo nulo en el borde

\[ m(K) \frac{ \rho_K^{n+1} - \rho_K^n }{\Delta t} = \sum_{ K|L \in \mathcal E_K \cap \mathcal E_{\mathrm{int}} } m(K|L) \frac{\rho_{L}^n - \rho_K^n}{ |x_L - x_K| } \]

También se pueden añadir otras condiciones: Dirichlet, Robin, …

Programando la ecuación de transporte

El paquete Triangulate.jl

Podríamos trabajar con cualquier tipo de mallas. Pero vamos a trabajar con mallas triangulares (para mantenerlas en Elementos Finitos).

Para la construcción de la malla vamos a utilizar el paquete Triangulate.jl. Es una interfaz del algoritmo Triangle de Jonathan Richard Shewchuk.

Podremos una condición de masa sobre los diferentes triángulos, para que la malla sea pequeña.

using Triangulate, Printf, Plots

function create_mesh()
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([-1.0 -1.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    
    maxarea=1e-2
    area=@sprintf("%.15f",maxarea)
    
    (triout,vorout) = triangulate("pa$(area)DQ", triin)
    return triout
end;
mesh = create_mesh();
## Structure of the output
# mesh.pointlist[i,j] i-th coordinate of the j-th point
# mesh.trianglelist[i,j] = k means that mesh.pointlist[k] is the i-th vertex of the j-th triangle
# segmentlist means that the segment from segmentlist[1,j] to segmentlist[2,j] is on the boundary
centroids = zeros(2,numberoftriangles(mesh))
for j in 1:numberoftriangles(mesh)
    centroids[:,j] = sum( [ mesh.pointlist[:, k] for k in mesh.trianglelist[:, j] ] )/3
end
function plot_triangle!(j)
    x = mesh.pointlist[1,1:end]
    y = mesh.pointlist[2,1:end]
    plt = plot!(
        [ x[mesh.trianglelist[1:3,j]]; x[mesh.trianglelist[1,j]]],
        [ y[mesh.trianglelist[1:3,j]]; y[mesh.trianglelist[1,j]]],
        color=:blue,
        label="",
        legend = :outertopright,
        )
end

# Plot all triangles
plot()
for L in 1:numberoftriangles(mesh)
    plot_triangle!(L)
    scatter!([centroids[1,L]],[centroids[2,L]],label="",color=:green)
end
plot!()

Área

Como solo vamos a trabajar con un malla triangular, podemos calcular las áreas con la fórmula de Heron

function area_of_triangle(vertices::Vector{Vector{Float64}})    
    # Formula de Heron para el area del triangulo de vértices vertices[1], vertices[2] y vertices[3]
    l1 = sqrt( (vertices[1][1]-vertices[2][1])^2 + (vertices[1][2] - vertices[2][2])^2 )
    l2 = sqrt( (vertices[2][1]-vertices[3][1])^2 + (vertices[2][2] - vertices[3][2])^2 )
    l3 = sqrt( (vertices[3][1]-vertices[1][1])^2 + (vertices[3][2] - vertices[1][2])^2 )
    s = ( l1 + l2 + l3 ) / 2
    return sqrt( s * (s-l1) * (s-l2) * (s-l3)  )
end

areas = ones(numberoftriangles(mesh))
for k = 1:numberoftriangles(mesh)
    index_of_vertices_of_T = mesh.trianglelist[:,k]
    vertices = [ mesh.pointlist[:,index_of_vertices_of_T[i]] for i in 1:3 ]
    areas[k] = area_of_triangle(vertices)
end

Para mallas generales se puede utilizar que \[ m(K) = \int_K 1 \, dx = \frac 1 d \int_K (\mathrm{div}(x)) dx = \frac 1 d \int_{\partial K} x \cdot n_K \approx \frac 1 d \sum_{\sigma \in \mathcal E_K} m(\sigma) x_{\sigma} \cdot n_{K,\sigma} \]

donde \(x_\sigma \in \sigma\) y usar que, en dimensión 2, \(m(\sigma)\) es sólo la longitud de un segmento.

Vecinos

Vamos a almacenar las aristas

struct Edge
    index_of_start_point::Int # Index in mesh.pointlist
    index_of_end_point::Int   # Index in mesh.pointlist
    mid_point::Vector 
    triangles::Tuple{Int,Int} # Indexes in mesh.trianglelist to the sides
    unit_normal_exterior:: Vector{Real} # Unit vector exterior to triangles[1]
    length::Real
end
function unit_normal_vector(K::Int,edge::Edge)::Vector
    if edge.triangles[1] == K
        unit_normal = edge.unit_normal_exterior
    elseif edge.triangles[2] == K
        unit_normal = - edge.unit_normal_exterior
    else
        error("$edge is not an edge of $K")
    end
end;

Extraemos información adicional de la malla

Calculamos los vecinos y las aristas interiores, es decir, los triángulos con \(\sigma = K|L \in \mathcal E\)

# Store the interior edges
interior_edges = Vector{Edge}()

# Given a triangle, which triangle share which edges 
# neighbours[K] is a vector, whose elements are tuples (L,index_of_K|L)
neighbours = [Vector{Tuple{Int64,Int64}}() for _ in 1:numberoftriangles(mesh) ];

Ahora tenemos que rellenar estos vectores con la información de la malla

Code 
for K in 1:numberoftriangles(mesh)
    for L in K+1:numberoftriangles(mesh)
        if sum([ point in mesh.trianglelist[:,L] for point in mesh.trianglelist[:,K] ]) == 2
            
            indexes_of_points_in_edge = intersect( mesh.trianglelist[:,K] ,  mesh.trianglelist[:,L] )
            index_third_point_in_K = setdiff(mesh.trianglelist[:,K],indexes_of_points_in_edge)[1]
            index_third_point_in_L = setdiff(mesh.trianglelist[:,L],indexes_of_points_in_edge)[1]

            # Recover the points
            P = mesh.pointlist[:,indexes_of_points_in_edge[1]]
            Q = mesh.pointlist[:,indexes_of_points_in_edge[2]]
            R = mesh.pointlist[:,index_third_point_in_K]
            S = mesh.pointlist[:,index_third_point_in_L]

            edge_tangent = P - Q
            edge_length = sqrt( sum(edge_tangent.^2 ))
            normal = [edge_tangent[2],-edge_tangent[1]] / edge_length

            if sum( (S-R).* normal ) <= 0
                normal .*= -1.0 
            end

            edge = Edge(indexes_of_points_in_edge[1],
                        indexes_of_points_in_edge[2],
                        (P+Q)/2,
                        (K,L),
                        normal,
                        edge_length )
            append!(interior_edges, [edge])
            current_edge = length(interior_edges)
            
            append!(neighbours[K],[(L,current_edge)])
            append!(neighbours[L],[(K,current_edge)])
        end
    end 
end

Calculamos el flujo de \(K\) a \(L\)

function velocity_times_exterior_normal_on_edge(K::Int,edge::Edge)        
    sum( velocity( edge.mid_point ) .* unit_normal_vector(K,edge) ) 
end
velocity_times_exterior_normal_on_edge (generic function with 1 method)
function up_wind_flux(K,L,edge,ρ) 
    v = velocity_times_exterior_normal_on_edge(K,edge)
    return ρ[K] * max(v,0) + ρ[L] * min(v,0)
end
up_wind_flux (generic function with 1 method)

Avanzamos en tiempo y representamos la superficie con PyPlot (GR) no representa bien

Code 
using PyPlot # GR does not plot this well
pyplot(fmt=:png)

Inicializamos en \(\rho_0 = 1\) y dejamos correr el tiempo

velocity(x) = -x
Δt = 1e-2
T  = 1

ρ = ones(numberoftriangles(mesh))
δ = zeros(numberoftriangles(mesh))
anim = @animate for t in 0:Δt:T
    global ρ
    for K in 1:numberoftriangles(mesh) 
        δ[K] = 1/areas[K]*sum([
                interior_edges[index_edge].length*
                    up_wind_flux(K,L,interior_edges[index_edge],ρ) 
                for (L,index_edge) in neighbours[K]
                ]) 
    end 
    ρ -= Δt * δ 

    surface(centroids[1,:],centroids[2,:],ρ,title="t=$t",colorbar = false, zlims = (0,10))
end every 1; fps=5
gif(anim)

Bibliografía

References

Eymard, R., Gallouët, T. & Herbin, R. (1997) Finite volume methods. In: Handbook of numerical analysis.