Traducción por Fran Navarro para CAChemE.org (CC BY)
FEniCS es una herramienta libre, gratuita y fácil de usar para la resolución de ecuaciones en derivadas parciales (EDPs) o en inglés, partial differential equations (PDEs). FEniCS es un proyecto de colaboración para el desarrollo de conceptos innovadores y herramientas para la computación científica automatizada, con un enfoque particular en solución automatizada de ecuaciones diferenciales por el método de elementos finitos.
La metodología y el software desarrollado como parte del proyecto FEniCS están documentados en una serie de artículos de científicos y un libro.
Un buen punto de partida para los nuevos usuarios es el Tutorial de FEniCS. Este documento recogerá la pimera parte del mismo y se centrará exclusivamente en los fundamentos de FEniCS mediante su interfaz Python, ya que este es el método más sencillo para explorar FEniCS para principiantes. El objetivo del tutorial (completo) es una ayuda para empezar con FEniCS a través de una serie de ejemplos sencillos que dan información sobre cómo:
Las matemáticas de los ejemplos se mantendrán simples para centrarse mejor en la funcionalidad de FEniCS y su sintaxis. Esto significa que principalmente se utilizará la ecuación de Poisson y la ecuación de difusión dependiente del tiempo como problemas modelo, a menudo con datos de entrada ajustados de tal manera que se obtiene una solución muy sencilla que se puede reproducir exactamente por cualquier método de elementos finitos estándar sobre una malla uniforme estructurada. Esta última propiedad simplifica en gran medida la verificación de las implementaciones. De vez en cuando se utilizarán ejemplos físicamente más relevantes para recordar al lector que el cambio de la EDP y condiciones de contorno a algo más real a menudo puede ser una tarea trivial.
FEniCS puede dar la sensación de requerir un conocimiento profundo de la matemática abstracta del método de los elementos finitos, así como familiaridad con el lenguaje de programación Python. Sin embargo, resulta que muchos lectores son capaces de adquirir los fundamentos de elementos finitos y de programación Python, ya que ir junto con este tutorial. Simplemente sigue leyendo y prueba los ejemplos. ¡Te sorprenderá lo fácil que es resolver EDPs con FEniCS!
Todos los ejemplos analizados en el siguiente están disponibles como archivos de código fuente de Python ejecutable en un árbol de directorios. El material de este tutorial es principalmente una traducción en español de la información recogida en su página web FEniCS en su versión 1.4.
Para realizar los ejercicios de este tutorial se requiere acceso a un equipo en el que esté instalado el software FEniCS. La sección Instalación FEniCS explica brevemente cómo instalar las herramientas necesarias. En Windows, FEniCS no se mantiene acutualizado a la última versión por lo que se recomienda instalar Ubuntu 12.04 LTS en una máquina virtual como VirtualBox.
El primer ejemplo se referirá al problema de Poisson: $$\begin{split}- \nabla^2 u(\pmb{x}) &= f(\pmb{x}),\quad \pmb{x}\mbox{ in } \Omega, \\ u(\pmb{x}) &= u_0(\pmb{x}),\quad \pmb{x}\mbox{ on } \partial \Omega\thinspace .\end{split}$$
Donde, $u(\pmb{x})$ es la función desconocida, $f(\pmb{x})$ es la función preescrita o término fuente, $\nabla^2$ es el operador de Laplace (tambien escrito en como $\Delta$), $\Omega$ es el dominio espacial, y $\partial\Omega$ es el contorno de $\Omega$. Una EDP estacionaria como esta, junto con un conjunto completo de las condiciones de frontera, constituye un problema de valores en la frontera, que debe establecerse con precisión antes de comenzar a resolver con FEniCS.
En dos dimensiones espaciales con coordenadas $x$ e $y$, podemos escribir la ecuación de Poisson como: $$- {\partial^2 u\over\partial x^2} - {\partial^2 u\over\partial y^2} = f(x,y)\thinspace .$$
La incógnita $u$ es ahora una función de dos vairiables $u(x,y)$, definida en un doominio bidimensional $\Omega$.
La ecuación de Poisson surge en numerosos contextos físicos, incluyendo la conducción térmica, electrostática, difusión de sustancias, la torsión de barras elásticas, el flujo de fluido no viscoso, y ondas. Por otra parte, la ecuación aparece en las estrategias de división numéricos de sistemas más complicados de las PDE, en particular, las ecuaciones de Navier-Stokes.
La solución de un problema físico con FENICS consta de los siguientes pasos:
Ahora vayamos más en detalle para los pasos 2-4 en detalle. La característica clave de FEniCS es que los pasos 3 y 4 resultan en código bastante sencillo, mientras que la mayoría de otros software para las EDPs requieren mucho más código y son técnicamente más difíciles de programar.
FEniCS hace que sea fácil resolver EDPs si los elementos finitos se utilizan para la discretización en el espacio y el problema se expresa como un problema variacional. Para los lectores que no están familiarizados con los problemas variacionales, se hará una breve introducción al tema pero se anima a leer un libro adecuado del método de elementos finitos.
El paso clave del procseo para convertir una EDP en un problema variacional es multiplicar la EDP por una función de $v$, integrar la ecuación resultante sobre $\Omega$, y llevar a cabo la integración por partes de los términos con derivadas de segundo orden. La función $v$ que multiplica la EDP en la literatura de los elementos finitos se llama una función de prueba (test function, en inglés). La función desconocida $u$ llamada función de aproximación (trial function, en inglés). La términos de test y funciones trial se utilizan en FEniCS también. Espacios funcionales adecuados se deben especificar para las funciones de test y de trial. Este tipo de espacios son bien conocidos para una EDP estándar que se plantean en la física y la mecánica.
Para este caso, lo primero que multiplicamos la ecuación de Poisson de la test function $v$ e integramos: $$-\int_\Omega (\nabla^2 u)v \, \mathrm{d}x = \int_\Omega fv \, \mathrm{d}x\thinspace .$$
Seguidamente integramos por partes el integrando con derivadas de segundo orden: $$-\int_\Omega (\nabla^2 u)v \, \mathrm{d}x = \int_\Omega\nabla u\cdot\nabla v \, \mathrm{d}x - \int_{\partial\Omega}{\partial u\over \partial n}v \, \mathrm{d}s ,$$
donde ${\partial u\over \partial n}$ es la derivada de $u$ en la dirección normal hacia el exterior en contorno. Se requiere que la test function $v$ desaparezca en las partes del contorno donde $u$ es conocido, que en el presente problema implica que $v = 0$ en todo el contorno $\partial\Omega$. El segundo término del lado derecho de la última ecuación, por tanto, desaparece. De esto se deduce que
$$\int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x = \int_{\Omega} fv \, \mathrm{d}x \quad .$$Esta ecuación se supone que debe mantener durante todo $v$ en algún espacio la función $\hat{V}$. La trial function $u$ radica en algun (posiblemente diferente) espacio funcional $V$. Se suele decir que la última ecuación es la forma débil (weak form, en inglés) del problema de contorno original, que consiste en la EDP $-\nabla^2u=f$ y la condición de contorno $u = u0$.
La declaración adecuada de nuestro problema variacional ahora sigue como: Encontrar $u \in V$ tal que
$$\int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x = \int_{\Omega} fv \, \mathrm{d}x \quad \forall v \in \hat{V}.$$Los espacios test y trial spaces $\hat{V}$ and $V$ se definen en el problema como: $$\begin{split}\hat{V} &= \{v \in H^1(\Omega) : v = 0 \mbox{ on } \partial\Omega\}, \\ V &= \{v \in H^1(\Omega) : v = u_0 \mbox{ on } \partial\Omega\}\thinspace .\end{split}$$
En resumen, $H^1 (\Omega)$ es el bien conocido en matemáticas espacio de Sobolev que contiene las funciones $v$ tal que $v^2$ y $||\nabla v||^2$ tienen integrales finitas $\Omega$. La solución de la EDP subyacente debe estar en función de un espacio donde también las derivadas son continuas, pero el espacio de Sobolev $H^1 (\Omega)$ permite funciones con derivadas discontinuas. Este requisito de continuidad más débil de $u$ en la forma variacional, provocada por la integración por partes, tiene grandes consecuencias prácticas cuando se trata de la construcción de los elementos finitos.
Para resolver la ecuación de Poisson numéricamente, necesitamos transformar el problema variacional continua a un problema variacional discreto. Esto se realiza mediante la introducción espacios (test y trial) de dimensión finita, a menudo denotadas como $\hat{V}_h\subset\hat{V}$ y $V_h\subset{V}$. El problema variacional discreto sige: Encuentra $u_h \in V_h \subset V$ tal que:
$$\int_{\Omega} \nabla u_h \cdot \nabla v \, \mathrm{d}x = \int_{\Omega} fv \, \mathrm{d}x \quad \forall v \in \hat{V}_h \subset \hat{V}\thinspace .$$La elección de $\hat{V}_h$ y $V_h$ se desprende directamente de la clase de elementos finitos queremos aplicar en nuestro problema. Por ejemplo, la elección del elemento triangular lineal conocida con tres nodos implica que $\hat{V}_h$ y $V_h$ son los espacios de todas las funciones lineales a trozos de una malla de triángulos, en donde las funciones en $\hat{V}_h$ son cero en el contorno y en los $V_h$ igual $u_0$ en el contorno.
En la literatura matemática de problemas variacionales se escribe $u_h$ para la solución del problema discreto y $u$ para la solución del problema continuo. Para obtener (casi) una relación de uno a uno entre la formulación matemática de un problema y el programa FEniCS correspondiente, vamos a utilizar $u$ para la solución del problema discreto y el $u_e$ para la solución exacta del problema continuo si queramos distinguir explícitamente entre las dos. En la mayoría de los casos, vamos a presentar el problema de la EDP con $u$ como desconocido, derivar una ecuación variacional a $(u, v) = L (v)$ con $u\in V$ y $u\in V$ y $v\in \hat V$, para luego simplemente discretizar el problema diciendo de elegir espacios de dimensión finita de $V$ y $\hat V$. Esta restricción de $V$ implica que $u$ se convierte en una función de elementos finitos discretos. En la práctica, esto significa que nos volvemos nuestro problema EDP en un problema variacional continuo,creamos una malla y especificamos un tipo de elemento, y luego dejamos $V$ que corresponda a esta malla y la elección elemento. Dependiendo de si $V$ es infinito o finito-dimensional, $u$ será la solución exacta o aproximada.
Es conveniente introducir la siguiente notación unificada para formas débiles lineales: $$a(u, v) = L(v)\thinspace .$$
En el problema presente tenemos que: $$\begin{split}a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v \, \mathrm{d}x, \\ L(v) &= \int_{\Omega} fv \, \mathrm{d}x\thinspace .\end{split}$$
En la literatura matemática, $a (u, v)$ que se conoce como una forma bilineal y $L (V)$ como una forma lineal. Tendremos que, para cada problema lineal, resolver e identificar los términos con la $u$ desconocida y recogerlas en $a(u, v)$, y de manera similar, recoger todos los términos con funciones sólo conocidos en L (v). Las fórmulas para $a$ y $L$ son luego porgramadas directamente en el en FEniCS.
En resumen, antes de hacer un programa FEniCS para resolver una EDP, se deben realizar dos pasos:
El problema a resolver hasta el momento tiene un dominio $\Omega$ general y funciones generales $u_0$ y $f$. Para nuestra primera aplicación debemos decidir sobre las opciones específicas de $\Omega$, $u_0$ y $f$. Será prudente escoger un problema específico en el que podemos comprobar fácilmente que la solución calculada es correcta. Comencemos con la especificación de una solución exacta: $$u_{\rm e}(x, y) = 1 +x^2 + 2y^2$$ en un dominio 2D. Mediante la inserción de la ecuación anterior en nuestro problema de Poisson, encontramos que la $u_e (x, y)$ es una solución si $$f(x,y) = -6,\quad u_0(x,y)=u_{\rm e}(x,y)=1 + x^2 + 2y^2,$$ independientemente de la forma del dominio. Elegimos aquí, para simplificar, el dominio que sea un cuadrado con longitud unidad $$\Omega = [0,1]\times [0,1] .$$
la razón para especificar la ecuacioón de la función exacta $u_e$ escogida es que el método de los elementos finitos, con un dominio rectangular dividido uniformemente en elementos triangulares lineales, reproducirá exactamente un polinomio de segundo orden en los vértices de las células, independientemente del tamaño de los elementos. Esta propiedad nos permite verificar la ejecución mediante la comparación de la solución computarizada, llamada $u$ en este documento (excepto cuando se configura el problema EDP), con la solución exacta, denotado por $u_e: u$ debe ser igual a $u$ en precisión de máuina en los nodos. Problemas de prueba con esta propiedad se construirán con frecuencia en el tutorial de la web de FEniCS.
Un programa FENICS para resolver la ecuación de Poisson en 2D con las opciones dadas de $\Omega$, $u_0$ y $f$ puede el que sigue:
In [ ]:
"""
FEniCS tutorial programa demo: Ecuación de Poisson con condiciones Dirichlet.
Ejemplo más simple de cálculo computacinal y visulación con FEniCS.
-Laplace(u) = f en un cuadrado de lado unidad.
u = u0 en el contorno.
u0 = u = 1 + x^2 + 2y^2, f = -6.
"""
# Los comentarios en Python se establecen con almohadillas
# esta información escrita se ignorará en la ejecuación
# Cargamos la librería de dolfin (FEniCS)
from dolfin import *
# Creamos la malla y definimos el espacio de fnciones
mesh = UnitSquareMesh(6, 4)
V = FunctionSpace(mesh, 'Lagrange', 1)
# Se definen las condiciones de contorno
u0 = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]')
def u0_boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u0, u0_boundary)
# Se define el problema variacional
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx
# Calcula la solución
u = Function(V)
solve(a == L, u, bc)
# Representa/dibuja la solución
plot(u)
# Representa/dibuja la malla
plot(mesh)
# Volcado de solución a presentar en formato VTK (opcional)
file = File('poisson.pvd')
file << u
# Mantiene la solución representada y la hace interactiva
interactive()
El código completo se puede encontrar en el d1_p2D.py archivo en el directorio de stationary/poisson. Para ejecutarlo basta con escribir en la línea de comandos
$ python run d1_p2D.pySe recomienda seguir la lectura del Tutorial de FEniCS así como un ejemplo de un problema real.