Matrices de coeficientes de los sistemas de OpenFOAM. Salida y representación visual con Python.

Ha pasado bastante tiempo desde la última vez que vine por aquí, pero ha sido una de esas épocas en las que enlazas unas cosas con otras y estás tan metido en un proyecto que no quieres dedicar un minuto a otra cosa.

Como ya os conté, lo que me quita ahora el sueño es el código de OpenFOAM, tanto crear solvers nuevos como conocerlo a fondo -no hace falta lo segundo para lo primero, creedme.

La idea ahora mismo es formatear los datos que crea OpenFOAM para poder usarlos en librerías de álgebra optimizada. Sinceramente, no se si esperar conseguir mejoras de tiempo o no. Lo hago porque se me ha metido en la cabeza y porque es una buena manera de soltarse, además el siguiente paso es exprimir CUDA en ello…

Bien, un par de puntos que hay que entender cuando digo que hay que formatear los datos de salida.


  • OpenFOAM genera unas matrices, mejor dicho unos objetos matrices, a partir de las clases básicas contenidas en su código. En la documentación podréis ver que hay scalarSquareMatrix, fvVectorMatrix, lduScalarMatrix… y un sinfín. No nos asustemos, el programita está creado usando templates para absolutamente todo, una forma bastante buena de ahorrar código…
  • Por lo general cada objeto no es solamente una matriz, sino que contiene al sistema completo de ecuaciones (clase fvMatrix) o la descomposición ldu de la matriz de coeficientes. Esto hace que debas tener un máster en punteros y en herencia de clases, pero es efectivo… =)
  • Para más inri, OpenFOAM usa un formato propio para almacenar matrices dispersas -que en dinámica de fluídos suelen ser todas- de esta manera sólo se almacenan los términos no-nulos reduciendo la memoria empleada. El truco para esto está en que otro array almacena las posiciones de los coeficientes no-nulos dentro de la matriz. De nuevo un máster para reconstruir la dichosa matriz.
  • Las matrices triangular inferior y superior tienen el mismo número de elementos no-nulos que caras tenemos en la malla. El número de elementos de la diagonal es igual al número de celdas. Esto hace que la matriz tenga una pinta interesante si pudiésemos verla de algún modo […]

Cuando lo anterior se comprende no es tan complejo obtener una salida con los coeficientes dispuestos en la forma correcta.

 

Una vez que tenía ese fichero me he hartado a hacer scroll para ver dónde estaban los términos no-nulos y que pinta tenía la matriz, así que he decidido representar la matriz acordándome de mis amigos de Pibonacci -gracias por enseñarme a cargar matrices desde un .dat de forma civilizada- y con un pequeño script en Python tenemos lo que sigue:

Es la forma más simple para ver si realmente la matriz tiene la pinta que debería. Algo a tener en cuenta es que la imagen está “invertida” respecto a los datos; la diagonal debe ser la bisectriz del segundo y cuarto cuadrante, no del primero y tercero, pero la representación es correcta.

Como la mayor parte de los datos son nulos no veo sentido a tirar de estadística para ver cosas, podeís comprobarlo empleando la matriz de correlaciones si queréis.

Os dejo el -corto y sencillo- script:

import numpy as np
from math import sqrt
from numpy.random import rand
from pylab import pcolor, show, colorbar, xticks, yticks

print 'Script para plotear una matriz contenida en un fichero de texto plano'

vDatos = np.loadtxt('dataMatrix.dat')

pcolor(vDatos)
colorbar()
show()

Dentro de nada os cuento algo más, sobre todo de BLAS y LAPACK. Siento no poder dejar aquí el código del solver.

Gracias a la Fundación Centro de Supercomputación de Castilla y León por los recursos.

Esta entrada fue publicada en Aerodinamica, C/C++, OpenFOAM, python. Guarda el enlace permanente.

3 respuestas a Matrices de coeficientes de los sistemas de OpenFOAM. Salida y representación visual con Python.

  1. Juanlu001 dijo:

    De nada!🙂 Y por cierto, puedes probar `imshow` para ahorrarte el rollo del eje invertido.http://matplotlib.org/examples/pylab_examples/pcolor_demo2.html

  2. pybonacci dijo:

    Por cierto (y perdona por el comentario doble), no sé qué formato de matriz dispersa usa OpenFOAM exactamente, pero SciPy tiene un módulo para matrices dispersas y aquí lo expliqué para el caso de matrices tridiagonales:https://pybonacci.wordpress.com/2012/05/10/como-crear-una-matriz-tridiagonal-en-python-con-numpy-y-scipy/Un saludo!

  3. No tiene importancia lo del comentario! Tienes razón, es buena idea lo de usar 'imshow' sabía que estaba por ahí pero como me hacía falta algo rápido he tirado de memoria, es un defecto! =)Sí, el módulo de matrices dispersas de SciPy pensaba usarlo para hacer cálculos rápidos, (realmente me siento mucho más cómodo escribiendo código que en C++ cuando hay temas de cálculos de por medio) pero no me da ninguna ventaja con OpenFOAM porque las matrices se guardan comprimidas y ya que hago el trabajo de "decodifciarlas" prefiero hacerlo directamente para C que es con lo que tengo que trabajar después.A la hora de medir tiempos quiero comparar distintos paquetes de álgebra y Numpy no se va a quedar fuera… ya os iré contando. Gracias por todo de nuevo!

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s