Lo que he aprendido: leer matrices MMF en Julia

Siguiendo con mi aprendizaje de Julia como sustituto a Matlab, he migrado una función que tenía para leer las matrices de masa y rigidez que exporto de Ansys para luego hacer movidas locas. Otro día os cuento el por qué de la transición (ejem, licencia, ejem) y cómo podéis uniros a mí en el lado oscuro.

Ale, al tema. Primero os voy a contar cómo exportar de Ansys en un par de párrafos, os doy un poco de contexto sobre escribir matrices en archivos de texto y luego vamos al grano.

Exportar matrices de Ansys

Para obtener las matrices de masa y rigidez de Ansys en primer lugar
exportamos el archivo input del modelo con: Tools > Write input file. Luego, a ese input le añadimos antes de la línea /wb,file,end este trocico:

! Gets Stiffness Matrix
! D for double precision real number, C for complex
*SMAT, MatK, D, import, full, file.full, stiff
*export, MatK, mmf, matkMMF.txt ! Exports Stiffness as MMF format

! Gets Mass Matrix
*SMAT, MatM, D, import, full, file.full, mass
*export, MatM, mmf, matmMMF.txt ! Exports Stiffness as MMF format

Como hemos usado SMAT exportará una matriz dispersa en formato MMF,
si quisiéramos obtener una matriz densa, tendríamos que usar DMAT.

Aparte, hacemos un Mechanical APDL, le metemos el input que hemos
modificado y calculamos. Así, nos creará los archivos matmMMF.txt y
matkMMF.txt en PROYECTO_files\dp0\APDL\ANSYS\

El formato MMF

De esta manera Ansys nos fabrica unos archivos de texto con un formato estándar para escribir matrices que en mi caso tiene esta pinta:

%%MatrixMarket matrix coordinate real symmetric 
%============================================================
% Creator: ANSYS Inc
% 
% This is a 24 x 24 Matrix
% There are 156 non-zero terms
%============================================================
24 24 156
1 1 9.541268758307953E+09
3 1 -3.747894585115073E+09
5 1 -2.807219846570191E+06
7 1 -5.498596390076558E+09
9 1 2.863561235766106E+08

Así sabemos que nos da la matriz en forma de coordenadas (posición, posición, valor), que es real y simétrica. Lo de que sea simétrica es clave saberlo, porque significa que solo va a escribir los valores del triángulo inferior, que no está por gastar espacio a lo tonto.

Otro tema es que la primera fila de números hace referencia al tamaño de la matriz y al número de términos que no son cero. Esto nos justifica el haber exportado como matriz dispersa, como matriz densa nos hubiera ocupado casi cuatro veces más.

Leer el archivo de texto y montar la matriz en Julia

A la hora de leer en Julia, he usado la función readdlm del paquete DelimitedFiles que me permite decirle qué separa los números, qué las líneas y cargarme los comentarios, en este caso las líneas que empiezan por un porcentaje. También he necesitado los paquetes SparseArrays y LinearAlgebra para poder fabricar la matriz dispersa y que Julia me cree una matriz simétrica a partir de la matriz triangular inferior.

Queda bastante finolis el código (¡y sin saber programar!):

using DelimitedFiles, SparseArrays, LinearAlgebra

function readMMF(file)

    data = readdlm(file,' ','\n'; comments=true, comment_char='%')

    N = data[1,1]; # number of elements

    # Data for sparse matrix
    i = data[2:end,1]
    j = data[2:end,2]
    v = data[2:end,3]

    A = sparse(i, j, v, N, N)
    # Create symmetric matrix from lower
    return Asym = Symmetric(A, :L)
end

Bonus: calcular valores y vectores propios

Con esto ya se pueden leer las matrices, pero a mí me mola el más difícil todavía y quería comprobar que los autovalores me dan igual que en Ansys (porque si no menuda liada). Después de leerme como media docu de Julia, porque están migrando paquetes y cosas, llegué a la conclusión de que tenía que usar el paquete Arpack para que me calculase los autovalores de mis supermatrices dispersas iterativamente.

Escribí un rollo así:

using Arpack

K = readMMF("K.txt")
M = readMMF("M.txt")

# Compute the first 15 eigenpairs
lambda, phi = eigs(K, M; nev=15, v0=zeros((0,)), which=:SR)
omega = sqrt.(filter(x -> x >= 0, lambda))

Aquí debo hacer unas cuantas notas:

  • Estoy calculando con la opción :SR porque me interesan los autovalores de menor parte real, es decir las primeras frecuencias naturales del sistema.

  • Le doy un vector inicial v0 porque si no empieza cada vez de uno aleatorio y si calculo 17 veces tengo 17 resultados parecidos pero diferentes. Empiezo con un vector cuyos elementos son cero, algo bastante típico.

  • Filtro los resultados antes de calcular las frecuencias naturales porque mi modelo está libre y aunque los primeros valores propios deberían dar cero, es posible que den números muy pequeños pero negativos y tendría que definir los autovalores como complejos cuando no lo tienen que ser para mi caso.

¡Y listo! Os dejo con la que probablemente sea la entrada más técnica del blog, espero que alguien la entienda 😀

Referencias

Matrix Market Exchange Formats

*Export Stiffness Matrix from Ansys *

Docu de Arpack

3 pensamientos en “Lo que he aprendido: leer matrices MMF en Julia

  1. Pingback: Compilación: los artículos del 2018 | Onda Hostil

¡Opina sin miedo! (Puedes usar Markdown)

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 )

Google+ photo

Estás comentando usando tu cuenta de Google+. 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 )

Conectando a %s