Source: https://pybonacci.org/2015/01/20/fenics-resolucion-de-ecuaciones-diferenciales-en-python/
Timestamp: 2020-07-03 20:27:49+00:00

Document:
FEniCS: Resolución de ecuaciones diferenciales en Python – Pybonacci
2015-01-20 2019-11-17
Etiquetas:condaecuaciones diferencialesEDPsfenicspython
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.
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:
Automatización de la discretización.
Automatización de la solución discreta.
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.
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.
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
conda activate fenics27 # Lo activamos
conda install fenics mkl --channel juanlu001 # Instalamos FEniCS y MKL
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á! 😀
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)
"""Ecuación de Poisson.
# Tenemos que especificar el número de elementos
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)
V = FunctionSpace(mesh, 'Lagrange', 1) # Polinomios de orden 1
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.
\( \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:
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.
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)
"""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.
u0 = Expression('1 + x[0] * x[0] + 2 * x[1] * x[1]')
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
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 😉
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
print max(abs(u.vector().array())) # Array de valores numéricos
Y para terminar, la gráfica correspondiente:
plot(u) interactive()
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:
https://twitter.com/pybonacci/status/553932429720055808
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!
WELLS, Garth; MARDAL, Kent-Andre; LOGG, Anders. Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book. Springer, 2012.
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 🙂
anteriorIntroducción a Machine Learning con Python (Parte 1)
siguiente¿Realmente es feo matplotlib?
19 pensamientos sobre “FEniCS: Resolución de ecuaciones diferenciales en Python”
CAChemEorg el 2015-01-23 a las 20:33
¡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.
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 😉
Juan Luis Cano el 2015-01-23 a las 22:10
¡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
Alberto el 2015-01-24 a las 20:01
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.
Juanlu001 el 2015-01-25 a las 17:06
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 😛
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.
Andrea el 2015-02-23 a las 20:57
Thank you for this fenics recipe!
I tried to install fenics on Ubuntu 14.10 with:
conda install fenics –channel juanlu001
Error: Could not find some dependencies for fenics: eigen3
Juan Luis Cano el 2015-02-23 a las 21:01
Hello Andrea! I saw your issue in GitHub, let us continue the conversation there. Thanks!
Anderson Lincol Condori Paytan el 2015-05-19 a las 06:28
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
Juan Luis Cano el 2015-05-30 a las 14:38
¡Hola Anderson! Tengo muy poca experiencia en mallados pero te diría que probaras gmsh.
Gerardo Loaiza el 2016-04-16 a las 02:44
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.
Juan Luis Cano el 2016-04-18 a las 17:33
Hola Gerardo, FEniCS no trae generador de geometrías así que tendrías que usar gmsh u otros. ¡Un saludo!
Gerardo el 2016-04-29 a las 04:48
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.
Juan Luis Cano el 2016-05-03 a las 11:45
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.
Gerardo el 2016-05-04 a las 20:29
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.
Gerardo el 2016-05-25 a las 21:53
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.
Gerardo el 2016-05-25 a las 21:56
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};
Anónimo el 2016-09-21 a las 19:57
¿Este método de instalación de Fenics también funciona para mac?
Mario el 2016-11-09 a las 18:31
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?
Juan Luis Cano el 2016-11-10 a las 00:18
A los que preguntan por Mac: para instalar FEniCS lo mejor es utilizar Docker. https://fenicsproject.org/download/#docker-images-all-platforms-and-versions-hspace-dockerimage
Muchas gracias a todos por vuestros comentarios, os animamos a utilizar http://es.stackoverflow.com/ para dudas concretas sobre código 🙂 ¡Saludos!

References: Resolución 
 resolución 
 resolución 
 resolución 
 Resolución 
 resolución