FEniCS: Resolución de ecuaciones diferenciales en Python

Introducción

En este artículo os voy a presentar el proyecto FEniCS, una colección de bibliotecas escritas con interfaz en Python para la resolución de ecuaciones diferenciales por el método de los elementos finitos. FEniCS proporciona un método muy inteligente para automatizar los procesos más laboriosos de la solución de este tipo de ecuaciones, de forma que podemos atacar cualquier problema que se nos ocurra. Las posibilidades son inmensas y la documentación es bastante exhaustiva, así que aquí veremos una introducción a cómo resolver ecuaciones en derivadas parciales con FEniCS.

Placa rectangular simplemente apoyada resuelta con FEniCS
Placa rectangular simplemente apoyada resuelta con FEniCS

Sin duda alguna el mejor recurso para aprender en profundidad a manejar las diferentes bibliotecas que conforman el proyecto es "The FEniCS book", disponible gratuitamente en Internet y en papel en la web de Springer. Muchas partes de este artículo están inspiradas o directamente traducidas de él. Está siendo mi biblia durante estas semanas y lo recomiendo al 100 %.

En esta entrada se han usado python 2.7.9, numpy 1.9.1 y fenics 1.5.0.

Este artículo tiene licencia GNU Free Documentation License, versión 1.3.

Alcance y objetivos

El proyecto FEniCS busca automatizar la resolución de ecuaciones diferenciales, utilizando el método de los elementos finitos. Este método, que data de mediados del siglo XX y que hoy por hoy es tremendamente importante en la resolución de problemas ingenieriles, se puede automatizar en tres pasos:

  1. Automatización de la discretización.
  2. Automatización de la solución discreta.
  3. Automatización del control del error.

FEniCS se creó en 2003 y se ocupa fundamentalmente de la automatización de la discretización. Una vez que se ha discretizado el problema y se han ensamblado los correspondientes sistemas lineales, su resolución (es decir, la segunda tarea) se deja a bibliotecas especializadas como PETSc, UMFPACK y otras.

Componentes del proyecto FEniCS (un poco desactualizado)
Componentes del proyecto FEniCS (un poco desactualizado)

FEniCS proporciona una interfaz en Python para todos estos paquetes, la mayoría escritos en Fortran o C++, y eso le da una sencillez de uso inigualable. Pero primero tenemos que ver cómo lo instalamos.

Instalación

Debido al gran número de dependencias de FEniCS, instalarlo puede resultar un desafío para los menos pacientes. Los desarrolladores proporcionan paquetes para Ubuntu (olvídate de Windows), y hay diversas soluciones para instalarlo en otras distribuciones.

Estas vacaciones invernales he puesto mi granito de arena y he creado unos paquetes y recetas conda para que se pueda instalar FEniCS junto a la distribución Anaconda. Si tienes suerte, es posible que puedas instalarlo en tres comandos:


$ conda create -n fenics27 python=2.7  # Creamos un entorno apropiado
$ source activate fenics27  # Lo activamos
(fenics27) $ conda install fenics mkl --channel juanlu001  # Instalamos FEniCS y MKL
(fenics27) $ python ~/.miniconda3/envs/fenics27/share/dolfin/demo/documented/poisson/python/demo_poisson.py  # ¡Cruza los dedos!

Los he compilado con la extensión MKL que proporciona Continuum, así que necesitarás una licencia académica gratuita.

Si por cualquier motivo la instalación falla (¡abre una incidencia!) o quieres modificar el proceso de compilación, no tienes más que bajarte las recetas y modificarlas a tu gusto. En el repositorio hay instrucciones sobre cómo utilizarlas. No me voy a detener mucho sobre esto porque sería material para otro artículo... ;)

En última instancia no te cuesta nada bajarte Linux Mint o Ubuntu, configurarlo en una máquina virtual como VirtualBox y empezar a trabajar. ¡Vamos allá! :D

La ecuación de Poisson

Antes de lanzarnos al código, como gente seria y ordenada que somos vamos a empezar con las matemáticas. Vamos a seguir los primeros apartados del tutorial de FEniCS, y vamos a resolver la ecuación de Poisson:

\(\displaystyle -\nabla^2 u = f \quad u \in \Omega.\)

De las condiciones de contorno hablaremos en seguida. Lo primero y más importante es definir el dominio, y en este caso resolveremos la ecuación en un cuadrado unitario \(\Omega = [0, 1] \times [0, 1] \subset \mathbb{R}^2\).

Vamos a ir escribiendo el programa paso a paso: después de importar todo FEniCS (sí, está hecho así en todas las demos, ¡no hay que alarmarse!), creamos el dominio correspondiente utilizando UnitSquareMesh:


# coding: utf-8
"""Ecuación de Poisson.

"""
from dolfin import *

# Tenemos que especificar el número de elementos
mesh = UnitSquareMesh(10, 10)

Espacios de funciones

Tenemos que establecer a qué espacios van a pertenecer nuestras funciones. Este es un paso muy delicado porque si no lo hacemos bien nuestros resultados serán absurdos, o lo que es peor, ligeramente incorrectos.

El tema tiene bastante enjundia, así que solo mencionaré que los espacios que se utilizan en el método de los elementos finitos se llaman espacios de Sobolev y que para este caso nos interesan funciones incluidas en \(H^1(\Omega)\). Ejemplos de elementos finitos tipo \(H^1(\Omega)\) son los clásicos elementos polinómicos de Lagrange.

Veamos cómo se traduce esto a código Python:


V = FunctionSpace(mesh, 'Lagrange', 1) # Polinomios de orden 1

u = TrialFunction(V)
v = TestFunction(V)

Hemos definido un espacio de funciones \(V = H^1(\Omega)\) sobre el dominio con FunctionSpace, y hemos creado ya nuestra función incógnita y nuestra función. ¡Seguimos!

Formulación variacional y condiciones de contorno

Lo siguiente que tenemos que hacer es escribir la ecuación en forma variacional, débil o integral. Para ello:

  • Multiplicaremos por una función test \(v\),
  • integraremos sobre todo el dominio \(\Omega\),
  • aplicaremos integración por partes para relajar las condiciones de derivabilidad sobre la función incógnita \(u\), y
  • impondremos las condiciones de contorno.

En nuestro caso nos quedaría:

\(\displaystyle -\int_\Omega (\nabla^2 u) v dx = \int_\Omega f v dx\)

Y aplicando integración por partes (que en dimensiones superiores a uno se conoce como primera identidad de Green),

\(\displaystyle -\int_\Omega (\nabla^2 u) v dx = \int_\Omega \nabla u \cdot \nabla v dx - \oint_{\partial\Omega} \frac{\partial u}{\partial n} v ds\)

Y ahora viene la parte importante en la que imponemos las condiciones de contorno, en la integral sobre \(\partial\Omega\). Hay dos opciones:

  1. Condiciones tipo Dirichlet: \(u = u_0\). En este caso hacemos desaparecer la última integrar imponiendo que \(v = 0\) en \(\partial\Omega\). En este caso las condiciones se denominan esenciales.
  2. Condiciones tipo Neumann: \(\frac{\partial u}{\partial n} = g\). En este caso la condición de contorno se traslada de manera natural a la formulación variacional (¡solo hay que sustituir!), y por esto se llaman condiciones naturales.

Esta distinción es importante porque se reflejará de manera clara en el código. Si quieres una explicación más profunda cualquier libro de elementos finitos o de cálculo variacional lo explicará mejor que yo.

Para este caso vamos a imponer condiciones de contorno Dirichlet (esenciales): \(u = 1 + x^2 + 2 y^2\) en \(\partial\Omega\) utilizando DirichletBC. Necesitamos trasladar esta expresión al programa (mediante la clase Expression y una función que defina en qué parte del contorno se aplican. El código será el siguiente:


def boundary(x, on_boundary):
    """El parámetro on_boundary es verdadero si el punto está
    sobre el contorno. Podríamos hacer otro tipo de comprobaciones,
    pero en este caso basta con devolver este mismo valor.

    """
    return on_boundary

u0 = Expression('1 + x[0] * x[0] + 2 * x[1] * x[1]')
bc = DirichletBC(V, u0, boundary)

En FEniCS la forma variacional se expresa con un lenguaje simbólico de la siguiente manera:

\(\displaystyle a(u, v) = L(v)\)

Nuestro término fuente será \(f(x, y) = -6\). Al haber aplicado condiciones de contorno esenciales, la última integral será nula. El código es:


f = Constant(-6.0)  # Término fuente
a = inner(nabla_grad(u), nabla_grad(v)) * dx  # Miembro izquierdo
L = f * v * dx  # Miembro derecho

Ya estamos a punto de sacar una bonita gráfica ;)

Resolución del sistema

La forma más directa de resolver el sistema es utilizando la función solve. Vamos a reutilizar la variable u para la solución, de esta forma:


u = Function(V)
solve(a == L, u, bc)

print max(abs(u.vector().array()))  # Array de valores numéricos

Y para terminar, la gráfica correspondiente:


plot(u)
interactive()

¡Y este es el resultado!

Ecuación de Poisson

Conclusiones

Alguno se ha podido marear con tanta ecuación, y es normal. En realidad, esto es solo la punta del iceberg: el tutorial básico de FEniCS continúa particionando la frontera para poner diferentes condiciones de contorno en cada parte, utilizando formulaciones mixtas, calculando el error cometido con respecto a la formulación exacta... Y aún queda una treintena de capítulos en el libro.

Durante los últimos días he estado publicando gists de problemas resueltos con FEniCS. Por ejemplo, la ecuación de Helmholtz:

Lo fundamental es tener pericia matemática y sobre todo muchas ganas de sacarle jugo a FEniCS. Espero que os haya gustado el artículo, no dudéis en ponerme vuestras dudas y sugerencias en los comentarios y, si tiene éxito, habrá segunda parte ;) ¡Un saludo!

Referencias

  • WELLS, Garth; MARDAL, Kent-Andre; LOGG, Anders. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer, 2012.

Agradecimientos

Desde aquí quiero agradecer a Pepe Cercós el haberme animado a empezar mi proyecto fin de carrera y a Fernando Varas por haberme devuelto la fe en los profesores de mi escuela con sus aclaraciones y dedicación. A hombros de gigantes :)

Juan Luis Cano

Estudiante de ingeniería aeronáutica y con pasión por la programación y el software libre. Obsesionado con mejorar los pequeños detalles y con ganas de cambiar el mundo. Divulgando Python en español a través de Pybonacci y la asociación Python España.

More Posts - Website

Follow Me:
TwitterLinkedIn

19 thoughts on “FEniCS: Resolución de ecuaciones diferenciales en Python

  1. ¡Espectacular entrada, Juanlu!

    Pepe Cercós me recomendó también FEniCS hace un tiempo y mi intención era retomarlo este cuatrimestre, esta entrada ayuda bastante a los más novatos con esta biblioteca (como yo). Pero creo que la mejor aportación es la de poder integrarlo con Anaconda (¡b r a v o!).

    He estado mirando y parece que siempre se debe de describir el sistema de ecuaciones de forma variacional, ¿no?.
    Lo digo porque hacer esto es “meterse en harina” y, para gente que empiece, quizás la cantidad de conceptos que se debe de manejar sea muy alta. Afortunadamente la documentación (¡y libro!) que acompañan a FEniCS permiten, supongo, aprender estos conceptos y ponerlos en práctica.

    A mi siempre me ha parecido muy interesante FiPy (funciona en Windows también), que permite la resolución también de PDEs pero mediante el método de volúmenes finitos.

    http://www.ctcms.nist.gov/fipy/

    Digo esto porque FiPy permite resolver PDEs definiendo sólo los coeficientes de la forma canónica y acoplar varias físicas de forma “sencilla”:
    http://www.ctcms.nist.gov/fipy/documentation/FAQ.html

    Por supuesto, FEniCS es algo diferente pero creo que también es interesante dejarlo por aquí escrito para futuros lectores de Pybonacci ;)

    1. ¡Muchísimas gracias por el comentario Fran!

      Prácticamente me he limitado a traducir el primer apartado del tutorial de FEniCS. El problema ha sido, como tú dices, que «me he tenido que meter en harina», y hasta el cuello: la parte matemática me ha costado semanas de estrujarme la cabeza y ha sido ahora cuando he empezado a ver la luz. La parte positiva es que he aprendido muchísimo acerca de conceptos que tenía muy mal aprendidos durante la carrera, y eso está siendo enormemente gratificante.

      No tengo ni idea de cómo funciona el método de los volúmenes finitos, pero seguramente para gente que esté habituada a hacer CFD FiPy les vendrá muy bien. Otra biblioteca parecida a FEniCS pero con un enfoque más práctico y orientado a mecánica de sólidos es SfePy:

      http://sfepy.org/doc-devel/index.html

      ¡Un saludo! :)

  2. Interesante entrada,como siempre. Y eso que el tema es lo suficientemente complejo como para no ser nada fácil escribir algo no demasiado extenso sobre el tema. La verdad es que la cantidad de conocimientos matemáticos que hay que tener para hincarle el diente al programa puede echar para atrás a muchos,pero las ecuaciones diferenciales son fundamentales en todos los ámbitos científico-técnicos y tarde o temprano toparemos con ellas. A mí personalmente me encantaría una entrada que explicase los distintos enfoques que hay para resolver EDP de forma numérica (incluso en qué casos hay una solución exacta y no sería necesario), con sus particularidades, ventajas y problemas, pero no sé si eso encajaría exactamente en la filosofía de Pybonacci, que prioriza (como no puede ser de otra manera) la generación de código por encima de temas mas teóricos.

    1. Muchas gracias por el comentario, Alberto. Desgraciadamente el tema que comentas, al margen de si entra o no dentro de lo que podríamos escribir en este blog, está bastante por encima de mi nivel de conocimientos :) Ojalá llegue algún día a ser capaz de materializar un artículo así. Mientras tanto, habrá más artículos de Python, que también se me da fatal pero por lo menos es más fácil :P

      ¡Un saludo!

  3. Sobre los programas para tratar EDP en Python si parece claro que los principales son FEniCS, FiPy y SfePy, como bien indicáis.
    Sin mas,felicidades por la entrada y por supuesto animarte a hacer no una, sino muchas mas sobre el tema, por su importancia y relevancia.
    ¡Un saludo!.

  4. Hi Juanlu,

    Thank you for this fenics recipe!

    I tried to install fenics on Ubuntu 14.10 with:
    conda install fenics –channel juanlu001
    and I get
    Error: Could not find some dependencies for fenics: eigen3

    Do you have any suggestion?

    Andrea

  5. Hola Juan Luis muy buen aporte y agradeciéndole por todos los aportes que nos brindas, y a la vez quisiera hacerles una consulta por favor si podrían orientarme.
    Cómo puedo obtener este tipo de mallado o resultados bidimensionales con python similares a las imágenes que adjunto:
    Las 3 primeras son desarrolladas por el método de volúmenes finitos y el ultimo enlace es una imagen de una red de agua todos ellos realizados en matlab. por favor espero puedan orientarme de como obtener estos tipos de gráficos para diferentes áreas irregulares, de antemano muchas gracias………….
    http://subefotos.com/ver/?4df8faa5042f30d8d625fe4148b8c1b5o.jpg
    http://subefotos.com/ver/?d47542b2fad35ca3858fcd6637a866eco.jpg#codigos
    http://subefotos.com/ver/?feda48fb7f3d053019cf9192274509c2o.jpg#codigos
    http://subefotos.com/ver/?11cee28e9f7a489009ca33a2146b63e9o.jpg#codigos

  6. Buenas noches. También me gusta fenics. He hecho cálculos con freefem++ y es bastante bueno. En fenics no he envontrado (como en freefem) la manera de hacer geometrias complejas. Por ejemplo un rectangulo con.un hueco circular por dentro. Si alguien conoce la forma sería genial que me indiquen en el fenicss book dónde leerlo.

    Gerardo desde Colombia.

  7. Gracias. GMSH es bastante bueno. El problema es que las mallas tienen extensión .mesh y fenics acepta .xml. He leido que la solución e susar dolfin-convert.py pero no logro que me funcione. Se supone que la sintaxis es:
    dolfin-convert miarchivo.mesh miarchivo.xml pero me sale un error: SyntaxError: invalid syntax. Si sabes cómo arreglar esto te lo agrdezco. Gerardo.

    1. Hola Gerardo, la verdad es que no entiendo a qué viene ese error y no encuentro referencias en Google, así que me temo que tendrás que preguntar a los desarrolladores. Un saludo.

      1. Buenas tardes. Bueno, probelam con dolfin-covert ya lo tengo resuelto. Por si a alguien le es útlil le eanexo mi solución (no es elegante, pero sirve): Dentro de FEniCS cree una subcarpeta (en mi caso: “C:\FEniCS\share\dolfin\demo\mis_programas”) y en ella ubiqué mi_archivo.py, mi_malla.msh, dolfin-covert.py y un .bat para ejecuatr cmd dentro de la carpeta. Una vez ejecutado el ,bat y linea d ecomandos escribo la orden que no me funcionaba:”dofin-convert.py mi_malla.msh mi_malla.xml” y es archivo .xml ya aparce y desde luego que se importa al scrpi d epython sin probela. La verdad no sé porqué no puedo ejecutar esa orden desde el IDL, pero ya no me importa. Me queda por averiguar si se pueden exportar a la malla las propiedades geométrica del malla hecha en GMSH. Feliz tarde.

  8. Buenas tardes. Aún tengo dificultades con dolfin-convert.py.
    He intentado n-veces importar la spropiedades físicas de una
    malla para utilizar condiciones de frontera con se explica en
    “http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/documented/bcs/python/documentation.html”.
    Seguí la sinstrucciones de los foros (renombrar “Physical groups”, etc. Pero nada.
    El archivo .geo es:<>
    EL archivo .msh es:<>

    Mallo, guardo la malla con circulos.msh “dolfin-conver.py circulos.msh circulos.xml” pero no me parece en el archivo .xml “cell_index” que según he leído es lo que me permite establecer las condicioens de frontera. Gracias.

  9. No cargó la URL: El archivo es:
    Point(1) = {0, 0, 0, 1}; Point(2) = {1, 0, 0, 1}; Point(3) = {-1, 0, 0, 1}; Point(4) = {0, 1, 0, 1}; Point(5) = {0, -1, 0, 1};
    Point(6) = {2, 0, 0, 0.5}; Point(7) = {-2, 0, 0, 0.5}; Point(8) = {0, 2, 0, 0.5}; Point(9) = {0, -2, 0, 0.5};
    Circle(1) = {8, 1, 6}; Circle(2) = {6, 1, 9}; Circle(3) = {9, 1, 7}; Circle(4) = {7, 1, 8}; Circle(5) = {4, 1, 2};
    Circle(6) = {2, 1, 5}; Circle(7) = {5, 1, 3}; Circle(8) = {3, 1, 4}; Line Loop(11) = {4, 1, 2, 3, -7, -6, -5, -8};
    Plane Surface(11) = {11}; Physical Line(1) = {1, 2, 3, 4}; Physical Line(2) = {5, 6, 7, 8};
    Physical Surface(0) = {11};

  10. Hola

    intente instalarlo en Anaconda con mac, en linux corre perfectamente, pero en el mac he intentado y no he tenido éxito.

    Sabes si hay que tener en cuenta algo más?

    gracias

Comments are closed.