MICSc : un código paralelo de dinámica de fluidos computacional basado en PETSc

MICSc: Un C´odigo Paralelo de Din´amica de Fluidos Computacional Basado en PETSc Proyecto Final de Carrera Autor: V´ıctor Gonz´alez Cort´es Dirigido por Jos´e E. Rom´an Escuela T´ecnica Superior de Ingenier´ıa Inform´atica, 1 de febrero de 2011 ii ´ Indice general Agradecimientos III 1. Introducci´on 1 1.1. Din´amica de fluidos computacional . . . . . . . . . . . . . . . 1 1.2. Paralelismo ............................ 2 1.3. Librer´ıas num´ericas . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4. Objetivos ............................. 5 2. Din´amica de fluidos computacional 7 2.1. Ecuaciones de Navier-Stokes . . . . . . . . . . . . . . . . . . . 7 2.2. Clasificaci´on de los fluidos . . . . . . . . . . . . . . . . . . . . 8 2.2.1. Subs´onico / Trans´onico / Supers´onico / Hipers´onico . 9 2.2.2. Compresible / Incompresible . . . . . . . . . . . . . . 9 2.2.3. Laminar / Turbulento . . . . . . . . . . . . . . . . . . 10 2.2.4. Estacionario / No estacionario . . . . . . . . . . . . . 10 2.2.5. Newtoniano / No newtoniano . . . . . . . . . . . . . . 11 2.3. Flujos turbulentos . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.1. DNS............................ 12 2.3.2. LES ............................ 12 2.3.3. RANS........................... 15 2.4. Sistema de ecuaciones lineales . . . . . . . . . . . . . . . . . . 15 2.4.1. Gradiente Conjugado . . . . . . . . . . . . . . . . . . 16 2.4.2. Gradiente Biconjugado . . . . . . . . . . . . . . . . . . 17 2.4.3. Residuo m´ınimo generalizado . . . . . . . . . . . . . . 19 2.5. Precondicionadores . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5.1. Jacobi........................... 20 2.5.2. ILU ............................ 21 2.5.3. Block Jacobi . . . . . . . . . . . . . . . . . . . . . . . 21 3. PETSc 23 3.1. Vectores.............................. 24 3.2. Matrices.............................. 26 3.3. Mallas estructuradas . . . . . . . . . . . . . . . . . . . . . . . 28 iii iv ´ INDICE GENERAL 3.4. Solvers lineales .......................... 29 3.5. Solvers nolineales ........................ 31 3.6. Integradores temporales . . . . . . . . . . . . . . . . . . . . . 31 4. MICSc 33 4.1. Discretizaci´on........................... 33 4.2. Entrada/Salida.......................... 35 4.2.1. Entrada.......................... 35 4.2.2. Salida ........................... 36 4.3. Sistema no lineal . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.3.1. Linealizaci´on . . . . . . . . . . . . . . . . . . . . . . . 37 4.3.2. M´etodo de Picard . . . . . . . . . . . . . . . . . . . . 38 4.3.3. Convergencia . . . . . . . . . . . . . . . . . . . . . . . 38 4.4. Estructuraci´on del c´odigo . . . . . . . . . . . . . . . . . . . . 39 4.4.1. Definici´on de las estructuras de datos . . . . . . . . . 39 4.4.2. Instancias y relaci´on de las estructuras de datos . . . . 46 4.4.3. Flujo de ejecuci´on . . . . . . . . . . . . . . . . . . . . 50 5. Desarrollos 54 5.1. Mantenimiento y optimizaci´on del c´odigo . . . . . . . . . . . 54 5.1.1. Estructuraci´on de librer´ıa y casos . . . . . . . . . . . . 54 5.1.2. Tratamiento din´amico de las variables . . . . . . . . . 57 5.1.3. Refactorizaci´on de los patches .............. 59 5.1.4. Refactorizaci´on de las propiedades . . . . . . . . . . . 62 5.2. Nuevas funcionalidades . . . . . . . . . . . . . . . . . . . . . . 63 5.2.1. Introducci´on de la viscosidad . . . . . . . . . . . . . . 63 5.2.2. C´alculo de medias . . . . . . . . . . . . . . . . . . . . 64 5.2.3. Modelo de van Driest . . . . . . . . . . . . . . . . . . 66 5.2.4. Documentaci´on . . . . . . . . . . . . . . . . . . . . . . 70 5.2.5. Salida de resultados . . . . . . . . . . . . . . . . . . . 71 5.2.6. Configure /Makefile ................... 74 6. Experimentos y resultados 77 6.1. Validaci´on LES: caso del canal peri´odico . . . . . . . . . . . . 77 6.2. Rendimiento del c´odigo paralelo . . . . . . . . . . . . . . . . . 79 6.2.1. Comparaci´on de solvers lineales y precondicionadores 80 6.2.2. Evaluaci´on de aceleraci´on y eficiencia paralela . . . . . 81 7. Conclusiones y trabajo futuro 83 7.1. Conclusiones ........................... 83 7.2. Trabajofuturo .......................... 84 A. Manual de usuario de MICSc 86 Agradecimientos Doy gracias a mi fam´ılia y amigos por el apoyo incondicional durante toda mi etapa de estudiante y porque sin ellos nada de esto habr´ıa sido posible. Tambi´en agradezco principalmente a Jos´e E. Rom´an y a Guillermo Palau que me hayan dado la oportunidad y todo el apoyo y material necesario para la realizaci´on del proyecto, as´ı como a Ana Cubero por permitirme formar parte de su trabajo y por su disposici´on y ayuda durante todo el trabajo. Adem´as, quiero dar las gracias, tanto al Vicerrectorado de Investigaci´on de la Universidad Polit´ecnica de Valencia y al Ministerio de Ciencia e Innovaci´on por promover y financiar el proyecto de investigaci´on RHELES en que se enmarca este proyecto final de carrera, como a todo el grupo GRyCAP por su apoyo, y en concreto a Andr´es Tom´as y Eloy Romero por su ayuda para este trabajo. v vi AGRADECIMIENTOS Cap´ıtulo 1 Introducci´on Este trabajo se sit´ua en el contexto del proyecto de investigaci´on RHELES, que tiene como principal objeto de an´alisis e investigaci´on la resoluci´on de sistemas de ecuaciones mediante algoritmos num´ericos y su aplicaci´on en el campo de la Din´amica de Fluidos Computacional (CFD por sus siglas en ingl´es, Computational Fluid Dynamics); concretamente, se centra en la resoluci´on iterativa de las ecuaciones de Navier-Stokes que rigen el comportamiento del fluido. En este cap´ıtulo se realiza una breve introducci´on a la CFD, as´ı como al paralelismo y las librer´ıas num´ericas, y finalmente se presentan los objetivos tanto del proyecto de investigaci´on RHELES como del proyecto final de carrera. 1.1. Din´amica de fluidos computacional La din´amica de fluidos computacional, consistente en la resoluci´on y el an´alisis de problemas relacionados con el flujo de fluidos mediante el uso de tecnolog´ıa computacional, tiene su comienzo en los a˜nos 60, cuando se reduc´ıa al ´ambito acad´emico y su visibilidad fuera de ´este era pr´acticamente nula. Sin embargo, en las ´ultimas d´ecadas ha experimentado una gran evoluci´on hasta convertirse en una de las herramientas industriales m´as importantes para el dise˜no y an´alisis de dispositivos de ingenier´ıa. Las causas principales de esta r´apida evoluci´on se deben, por un lado, al aumento de la potencia y el abaratamiento de los recursos computacionales; y, por otra parte, al desarrollo de algoritmos num´ericos para resolver de manera eficiente los problemas que se plantean en el contexto de la ingenier´ıa. De esta manera, se pueden destacar como hitos importantes en el desarrollo de la CFD, la publicaci´on en 1972 del primer m´etodo segregado para la resoluci´on de ecuaciones de flujo incompresible [1], o la aparici´on de aplicaciones comerciales como FLUENT y CFX entre otros. Posteriormente, en la ´ultima d´ecada, la CFD est´a experimentando un 1 2CAP´ ITULO 1. INTRODUCCI ´ ON nuevo desarrollo, en esta ocasi´on, propiciado por el aumento de la capacidad computacional y el progreso en los m´etodos iterativos para resolver sistemas de ecuaciones lineales. Adem´as, la mayor complejidad en el patr´on de llenado de la matriz de coeficientes, as´ı como el aumento de la exigencia debido al mayor uso industrial de la CFD, ha propiciado un replanteamiento de los algoritmos tradicionales. De esta manera, se requieren c´odigos eficientes a la par que robustos, capaces de resolver grandes problemas num´ericos en el menor tiempo posible, y minimizando el n´umero de par´ametros libres introducidos por el usuario. 1.2. Paralelismo El paralelismo o computaci´on paralela se define como la ejecuci´on concurrente de una misma tarea mediante diferentes unidades de proceso. Esta concurrencia se puede dar en distintos escenarios: dentro del propio procesador, en un ´unico computador o incluso en Internet. En el propio procesador existe la posibilidad de incluir varias unidades de c´alculo de manera que es posible ejecutar de manera simult´anea diversas operaciones b´asicas. Cuando la computaci´on paralela se da en Internet (o en cualquier otra red de interconexi´on local o global), aparece el concepto de tecnolog´ıa grid, que consiste en la uni´on de computadores en diferentes puntos del planeta conectados a trav´es de una red de forma que se comportan como un ´unico elemento de c´alculo. Por ´ultimo, el caso en que el paralelismo se produce en un mismo computador es precisamente el que se da a lo largo de todo este trabajo y presenta diferentes variantes. El paralelismo a nivel de computador consiste en la integraci´on de diferentes procesadores en un ´unico sistema y existen principalmente dos paradigmas. En primer lugar se da el paradigma de memoria compartida, en el que todos los procesadores tienen acceso a la misma memoria y por otro lado existe el paradigma de memoria distribuida en el que cada procesador tiene acceso exclusivo a una zona de memoria propia. En este segundo paradigma, que es el que se utiliza de manera exclusiva en este proyecto, se utiliza la interfaz de paso de mensajes (MPI por sus siglas en ingl´es, Message Passing Interface [2]). Este est´andar define las funciones necesarias para la comunicaci´on entre los distintos procesadores, ya que siendo la memoria distribuida se hace necesario alg´un mecanismo para la transmisi´on de datos entre diferentes procesadores. Las funciones principales que define MPI se pueden clasificar en tres grupos: Inicializar, administrar y finalizar comunicaciones: •MPI Init: Inicia una sesi´on MPI. Se debe llama siempre antes de cualquier otra funci´on de MPI. 1.2. PARALELISMO 3 •MPI Finalize: Termina una sesi´on MPI. Debe ser la ´ultima llamada a MPI del programa. •MPI Comm size: Obtiene el n´umero total de procesos. •MPI Comm rank: Obtiene el identificador (rank) del proceso. Comunicaciones punto a punto: •MPI Send: Env´ıa un dato a otro proceso. •MPI Recv: Recibe un dato de otro proceso. Para que se produzca la comunicaci´on ambos procesos deben llamar a las respectivas funciones, d´andose un bloqueo si alguna de las dos no se produce. Comunicaciones colectivas: •MPI Bcast: Difunde un dato de un proceso a todos los demas. •MPI Scatter: Distribuye desde un proceso pun conjunto de n datos entre los nprocesadores (un dato por procesador). Equivale a realizar desde p, por cada proceso, una llamada a MPI Send con el dato y destino correspondientes. •MPI Gather: Recibe en el procesador pun conjunto de ndatos desde los nprocesadores (un dato por procesador). Es la operaci´on inversa a MPI Scatter y equivale a realizar en cada proceso una llamada a MPI Send con el dato correspondiente y con destino el procesador p. •MPI Allgather: Recibe en todos y cada uno de los procesadores, ndatos desde los nprocesadores (un dato por procesador). Es equivalente a realizar una llamada a MPI Gather en cada uno de los procesadores. Por ´ultimo, tambi´en es interesante conocer las magnitudes que nos permiten medir el comportamiento paralelo de una cierta aplicaci´on. Es por esto que definiremos el concepto de speedup o aceleraci´on y eficiencia. La aceleraci´on se define como el factor de ganancia de un programa paralelo con respecto a la versi´on secuencial: Sp=t1 tp , donde t1es el tiempo del programa secuencial y tpes el tiempo del programa paralelo. Es com´un no disponer de la versi´on secuencial del algoritmo, por lo que generalmente se realiza la comparaci´on con respecto al programa paralelo ejecutado con un solo procesador. Tambi´en es importante destacar que el valor de esta magnitud viene condicionado por el n´umero de procesos y se pueden dar las siguientes situaciones: 10 CAP´ ITULO 2. DIN ´ AMICA DE FLUIDOS COMPUTACIONAL de la masa sin ning´un t´ermino que se corresponda con la presi´on de manera directa o indirecta. Puesto que generalmente esta ecuaci´on se utiliza para resolver la presi´on, esto puede llevarnos, tarde o temprano, a la aparici´on de un 0 en la diagonal del sistema de ecuaciones que debe ser resuelto de alguna manera; este problema se trata con mayor nivel de detalle en 4.1. 2.2.3. Laminar / Turbulento Se dice que un fluido es turbulento cuando est´a caracterizado por recirculaciones o “remolinos” que se mueven de manera ca´otica o aparentemente aleatoria. Se introduce en este contexto el concepto de Eddy, con el que nos referimos a estos torbellinos turbulentos que se crean en un fluido, por ejemplo, tras pasar junto a un obst´aculo. Debido a la complejidad del fen´omeno en s´ı y de que el proyecto da gran importancia a los flujos turbulentos, en 2.3 se comenta m´as detalladamente esta situaci´on. Por otro lado, un fluido en el que no se den estas circunstancias, donde todas las part´ıculas sigan una trayectoria relativamente uniforme, provocada exclusivamente por las fuerzas externas causantes del movimiento, se denomina laminar. Para la caracterizaci´on de un flujo turbulento o laminar es de gran utilidad el n´umero de Reynolds. Este n´umero es una magnitud adimensional que relaciona las fuerzas inerciales con las fuerzas debidas a la viscosidad de la siguiente manera: Re =ρV L µ donde ρes la densidad, Vla velocidad media del fluido, Lla longitud caracter´ıstica del sistema (por ejemplo, en el caso del flujo en una tuber´ıa, el di´ametro de dicha tuber´ıa) y µla viscosidad din´amica del fluido. De esta manera, conociendo el n´umero de Reynolds de un fluido, es posible aproximar el nivel de turbulencia del mismo. Sirva como ejemplo el caso del flujo a trav´es de una tuber´ıa, donde el flujo ser´ıa laminar cuando Re <2300, turbulento cuando Re >4000 y un estado de transici´on, donde ambas situaciones son posibles y dependen de otros factores, en el intervalo intermedio. Por ´ultimo, cabe destacar que, a pesar de que el fen´omeno de la turbulencia no simplifica ni modifica de manera directa las ecuaciones del fluido, s´ı que es de capital importancia su consideraci´on de cara a adoptar una de las distintas estrategias disponibles en los flujos turbulentos, tal y como se comenta m´as adelante. 2.2.4. Estacionario / No estacionario Los fluidos tambi´en pueden ser clasificados en funci´on de su comportamiento a lo largo del tiempo. De esta manera, un flujo cuyas propiedades no 2.2. CLASIFICACI ´ ON DE LOS FLUIDOS 11 var´ıen a lo largo del tiempo es estacionario, mientra que si estas propiedades cambian en el tiempo, se dice que el flujo est´a en r´egimen transitorio o que es no estacionario. Los fluidos turbulentos son no estacionarios por definici´on, ya que el propio fen´omeno turbulento propicia un movimiento continuo a trav´es de los Eddies o torbellinos. Sin embargo, Pope afirma que un determinado campo aleatorio U(x,t) es estad´ısticamente estacionario si todas las estad´ısticas permanecen invariables bajo un determinado paso de tiempo; esto es, todas las propiedades permanecen constantes en el tiempo [6]. De esta manera, tomando la media de los campos de inter´es en un flujo turbulento con un determinado paso de tiempo puede llevarnos a la invariabilidad de los mismos y, por tanto, a que el flujo pueda alcanzar un r´egimen estad´ısticamente estacionario. Por ´ultimo, la principal ventaja que plantean un flujo estacionario con respecto de un flujo no estacionarios es la eliminaci´on de todas las derivadas temporales de las ecuaciones de conservaci´on. As´ı pues, al margen de la ventaja debida a la estabilidad del problema para un flujo estacionario, el sistema resultante constar´a de una dimensi´on menos: el tiempo. 2.2.5. Newtoniano / No newtoniano Un fluido es newtoniano cuando la curva que relaciona el estr´es o fuerza aplicada al fluido y su deformaci´on (que se pone de manifiesto mediante un cambio en la velocidad del fluido) es una recta que pasa por el origen; es decir, est´an relacionados linealmente mediante una constante de proporcionalidad. Dicha constante es la viscosidad del fluido y se define de la siguiente manera: τ=µdu dy donde τes el esfuerzo cortante (o fuerza aplicada sobre el fluido de manera tangencial), µes la constante de proporcionalidad o viscosidad ydu dy es el gradiente de velocidad, perpendicular a la direcci´on del esfuerzo cortante. As´ı pues, un fluido cuya viscosidad sea constante ser´a newtoniano (agua, aire, vino, ...), mientras que en caso contrario ser´a no newtoniano (miel, pegamento, ...). Esto tiene una consecuencia directa en la parte derecha de las ecuaciones de conservaci´on del momento, puesto que parte de las fuerzas que act´uan sobre el fluido se pueden escribir en funci´on del esfuerzo cortante, que a su vez depende de una viscosidad constante. Como explican Ferziger y Peri´c [7], muy diferente es el caso de fluidos no newtonianos, donde la relaci´on entre el esfuerzo cortante y la velocidad se define mediante un conjunto de ecuaciones en derivadas parciales que aumentan en gran medida la complejidad del problema. Sin embargo, la mayor parte de los fluidos reales, y m´as concretamente los tratados en este proyecto, se comportan como fluidos newtonianos. 12 CAP´ ITULO 2. DIN ´ AMICA DE FLUIDOS COMPUTACIONAL 2.3. Flujos turbulentos Como se ha explicado en el apartado anterior, los flujos turbulentos se caracterizan principalmente por la aparici´on de recirculaciones de manera ca´otica o aparentemente aleatoria. As´ı, es necesario poder entender y predecir este comportamiento con la mayor exactitud posible con el fin de obtener buenos dise˜nos en la ingenier´ıa. Por otro lado, con el aumento en los requisitos de las aplicaciones industriales, se requieren cada vez mayores niveles de detalle y precisi´on, por lo que los m´etodos num´ericos se han vuelto imprescindibles en el estudio de flujos turbulentos. En las siguientes subsecciones definiremos con algo m´as de detalle en qu´e consisten las principales aproximaciones disponibles para la predicci´on de flujos turbulentos. Cabe destacar que, al margen de las que aqu´ı se exponen, existen m´as aproximaciones [8] pero que carecen de importancia para este trabajo y a la vez, exceden el alcance del mismo. 2.3.1. DNS El primer m´etodo es la simulaci´on num´erica directa (DNS por sus siglas en ingl´es, Direct Numerical Simulation). Se trata de la t´ecnica de simulaci´on de flujos turbulentos m´as precisa, ya que no se promedia ni se aproxima ning´un resultado, sino que se resuelven las ecuaciones de Navier-Stokes para todos los Eddies, independientemente de su tama˜no o escala. Constituye la aproximaci´on m´as simple desde el punto de vista conceptual y en cierto modo se puede considerar equivalente a la realizaci´on de un experimento en laboratorio con un flujo y dominio de las mismas caracter´ısticas. Sin embargo, para capturar todas las estructuras del fen´omeno turbulento, incluyendo las de menor escala, es necesario disponer de una malla con una densidad suficiente. Generalmente, esto llega a ser prohibitivo para la posterior ejecuci´on de la simulaci´on. No obstante, existen simulaciones bien conocidas, realizadas con la t´ecnica DNS, que suelen utilizarse para la validaci´on de resultados obtenidos con otra t´ecnica (u otro c´odigo) para el mismo problema. 2.3.2. LES El m´etodo de simulaci´on de grandes escalas (LES, por sus siglas en ingl´es Large Eddy Simulation) consiste en resolver o simular los Eddies de gran escala, mientras que los menores se modelan, pudiendo elegir entre diferentes modelos para este prop´osito. La idea b´asica tras esta t´ecnica es que los Eddies de mayor escala son mucho m´as energ´eticos que los m´as peque˜nos y, por tanto, contribuyen en mayor medida al transporte de las propiedades del flujo, tal y como explica Kolmogorov en su teor´ıa de 1941. De esta manera, sin entrar en detalles formales, el proceso consiste en dividir los campos 2.3. FLUJOS TURBULENTOS 13 del fluido (principalmente velocidades y presi´on) en la parte resuelta para las grandes escalas y la parte que debe ser modelada. As´ı, se resolver´an las ecuaciones de Navier-Stokes para las grandes escalas, pero se a˜nadir´a un t´ermino a estas ecuaciones que incluya la parte modelada. Es importante destacar que el filtro que determina qu´e se debe resolver y qu´e se debe modelar es, generalmente, la propia malla, por lo que los modelos se conocen como modelos SGS por sus siglas en ingl´es (SubGrid Scale). De esta manera, el t´ermino adicional que se introduce en las ecuaciones, y que representa la diferencia entre el campo total y el campo considerando ´unicamente las fen´omenos de mayor escala, es 1 ρ ∂τij ∂xj donde normalmente se utiliza la siguiente igualdad para el t´ermino τij: τij −1 3τkkδij = 2µtSij siendo Sij la velocidad de deformaci´on Sij =1 2∂ui ∂xj +∂uj ∂xi yµtla viscosidad turbulenta a simular. En los siguientes apartados trataremos el modelo de Smagorinsky para esta viscosidad turbulenta, y los modelos de van Driest y (de manera muy superficial) din´amico, que se basan en este primero. Modelo de Smagorinsky El modelo de Smagorinsky define la viscosidad turbulenta como µt=C2 Sρ∆2|S| donde ∆ es la ra´ız c´ubica del volumen finito, Ses p2SijSij yCSes la constante de Smagorinsky, propia del modelo y que suele tomar valores en el rango [0.065,1.1] dependiendo del tipo de flujo turbulento. Este modelo ha demostrado ser satisfactorio en la pr´actica. Sin embargo, plantea el principal problema de que la constante de Smagorinsky mantiene el mismo valor en todo el dominio, por lo que la viscosidad turbulenta no decrece conforme se considera el flujo m´as cercano a las superficies, donde el flujo tiende a ser laminar. Con el fin de atajar este inconveniente se crearon los siguientes modelos. 14 CAP´ ITULO 2. DIN ´ AMICA DE FLUIDOS COMPUTACIONAL Modelo de van Driest Como se ha comentado en el apartado anterior, la simple aplicaci´on del modelo de Smagorinsky cerca de las superficies nos lleva a una simulaci´on que se aleja de los resultados correctos a medida que la distancia a la superficie disminuye. Para contrarrestar esto se define la amortiguaci´on de van Driest [9], que consiste principalmente en multiplicar la constante de Smagorinsky (CS) por un factor de escala que disminuye dicha constante en funci´on de la distancia a las superficies de la siguiente manera: CS=CS(y) = CS·1−e(−y∗/A).(2.4) En esta expresi´on, Aes la constante de van Driest y suele tomar el valor 25, mientras que y∗(tambi´en conocido como n´umero de Reynolds turbulento), se calcula mediante y∗=d νrτw ρ,(2.5) donde des la distancia m´ınima a la pared, νes la viscosidad cinem´atica del fluido y τwel esfuerzo cortante. Adem´as, la viscosidad cinem´atica cumple que ν=µ ρ,(2.6) siendo µla viscosidad cinem´atica (par´ametro de entrada). Por ´ultimo, sabiendo que es posible calcular el esfuerzo cortante mediante τw=µvtan d,(2.7) donde vtan es la componente de la velocidad tangencial a la pared, es posible reescribir la ecuaci´on 2.5 como y∗=d νrτw ρ=d νrµ vtan d ρ =d νrν vtan d=rd vtan ν=sd ρ vtan µ(2.8) y sustituyendo en la ecuaci´on inicial (2.4), se obtiene CS=CS(d) = CS 1−e −rdρvtan µA2!.(2.9) Modelo din´amico Por ´ultimo, el modelo din´amico se dirige al mismo fin que el modelo de van Driest y plantea una soluci´on similar. En este caso, la constante de Smagorinsky pasa de ser constante en todo el dominio a ser un valor dependiente del espacio y el tiempo. As´ı, en cada instante de tiempo y cada 2.4. SISTEMA DE ECUACIONES LINEALES 15 volumen finito, el valor podr´a ser distinto. A pesar de que se considera un desarrollo necesario a corto plazo para el proyecto de investigaci´on, para el desarrollo de este trabajo no se ha necesitado este modelo, por lo que no se detallar´a formalmente de qu´e manera es posible calcular la constante de Smagorinsky en cada punto e instante de tiempo. 2.3.3. RANS Este m´etodo se basa en un conjunto de ecuaciones obtenidas promediando las ecuaciones de conservaci´on en el tiempo, que se conocen como ecuaciones Reynolds-Averaged Navier-Stokes (RANS). De esta manera, toda la componente no turbulenta es promediada; es decir, es considerada como parte de la turbulencia, lo cual introduce en las ecuaciones promediadas ciertos t´erminos que deben ser modelados. Adem´as es importante saber que, debido a la complejidad de los flujos turbulentos, es complicado que exista un modelo RANS capaz de representar todos los flujos, por lo que estos modelos y resultados deben considerarse como aproximaciones con aplicaci´on en la ingenier´ıa y no como leyes cient´ıficas. Por ´ultimo cabe destacar que, al margen de que est´a fuera del objetivo de este documento, esta t´ecnica no ha sido utilizada a lo largo de este trabajo, por lo que no entraremos a detallar los diferentes modelos. 2.4. Sistema de ecuaciones lineales Independientemente de la estrategia seguida para afrontar los flujos turbulentos (DNS, LES, RANS), la linealizaci´on del sistema de ecuaciones (Newton 3.5, Picard 4.3.1, ...) o la discretizaci´on (4.1), finalmente hemos de afrontar la resoluci´on de un gran sistema de ecuaciones lineales con una matriz de coeficientes de gran dispersi´on de la forma Ax =b. Para resolver este tipo de problemas existen dos aproximaciones fundamentales: los m´etodos directos y los m´etodos iterativos. Entre los m´etodos directos cabe destacar la factorizaci´on LU o eliminaci´on gaussiana, que es uno de los m´etodos m´as intuitivos y simples para resolver sistemas de ecuaciones lineales (ver [10] para m´as detalles). Sin embargo, este m´etodo produce llenado; es decir, a medida que el m´etodo avanza aumenta el n´umero de elementos no nulos de la matriz, por lo que en gran parte se pierde la ventaja de la dispersi´on de la matriz. Es por esto que el coste de este algoritmo resulta prohibitivo para la dimensi´on de nuestro problema y, en general, para problemas resultantes de discretizar ecuaciones de derivadas parciales. Por otro lado, se encuentran los m´etodos iterativos que se basan en dar una primera aproximaci´on arbitraria de la soluci´on e ir aproxim´andola de manera iterativa a la soluci´on exacta hasta que se cumple un cierto criterio de convergencia. Tal vez, la implementaci´on m´as simple de esta idea sea 16 CAP´ ITULO 2. DIN ´ AMICA DE FLUIDOS COMPUTACIONAL reescribir el sistema Ax =bcomo una iteraci´on del punto fijo. Esta iteraci´on se define sobre una funci´on f:R→Rde la siguiente manera: xn+1 =f(xn), n = 0,1,2, ... Esta funci´on genera una sucesi´on x0,x1,x2, ... que se espera que converja a un determinado valor x. Si esta sucesi´on es convergente es posible demostrar que esa xes un punto fijo de la funci´on; es decir, x=f(x). De esta manera, podemos reescribir el sistema Ax =bcomo x= (I−A)x+b. y definir la iteraci´on de Richardson de la siguiente manera: xk+1 = (I−A)xk+b. En general, nos referimos a cualquier m´etodo donde xk+1 =Mxk+c como m´etodo estacionario, ya que el paso de xkaxk+1 no depende de nada m´as que de el resultado anterior y de la matriz Momatriz de iteraci´on, que permanece invariable durante toda la resoluci´on. Sin embargo, en este trabajo no se hace uso de los m´etodos iterativos estacionarios, sino de aquellos basados en los subespacios de Krylov. Estos m´etodos no se basan en la existencia de una matriz de iteraci´on, sino en el subespacio generado por una matriz A de dimensi´on n×ny un vector b de tama˜no n. Este subespacio se genera mediante las primeras k(k < n) potencias de Aaplicadas sobre r0, esto es: Kk(A, r0) = span{r0, Ar0, A2r0, . . . , Ak−1r0} Partiendo de una soluci´on inicial arbitraria (generalmente, x0= 0), y haciendo uso de los subespacios de Krylov, los m´etodos que se exponen a continuaci´on, para una determinada iteraci´on k, tratan de minimizar alguna medida del error de la aproximaci´on calculada sobre el subespacio x0+Kk. 2.4.1. Gradiente Conjugado El m´etodo del gradiente conjugado (en adelante CG por sus siglas en ingl´es, Conjugate Gradient), se basa en la minimizaci´on de la siguiente funci´on cuadr´atica: φ(x) = 1 2xTAx −xTb. Es posible demostrar que esta funci´on alcanza un m´ınimo cuando Ax =b. De esta manera, en cada iteraci´on, el m´etodo realiza una b´usqueda unidimensional a lo largo de una determinada direcci´on skde manera que xk+1 =xk+αsk, donde αes el par´ametro que se debe determinar de forma 2.4. SISTEMA DE ECUACIONES LINEALES 17 que la funci´on φ(xk+αsk) se minimice. Este m´ınimo se produce cuando el nuevo residuo es ortogonal a la direcci´on de b´usqueda, es decir: rT k+1sk= 0. Llegados a este punto, es posible expresar el nuevo residuo en funci´on del par´ametro αde la siguiente manera: rk+1 =b−Axk+1 =b−A(xk+αsk) = b−Axk−αAsk=rk−αAsk y, por tanto, expresar αcomo: α=rT ksk sT k(Ask). Habiendo determinado la funci´on a minimizar, la actualizaci´on del vector soluci´on y el c´alculo del par´ametro a lo largo de la direcci´on de b´usqueda, s´olo queda determinar c´omo obtener la nueva direcci´on de b´usqueda. Desde el punto de vista conceptual, esta nueva direcci´on deber´ıa ser ortogonal a todas las direcciones anteriores, haciendo uso de la ortogonalizaci´on de Gram-Schmidt o similar. Sin embargo, esta ortogonalizaci´on requerir´ıa mayor almacenamiento para todas las direcciones utilizadas, as´ı como un coste computacional prohibitivo. Es por esto que en la pr´actica se opta por que la nueva direcci´on sea A-ortogonal (o conjugada, de ah´ı el nombre del m´etodo) con respecto a la direcci´on anterior. De esta manera, el coste computacional y de almacenamiento es m´ınimo y el algoritmo resultante es: Algoritmo 1 Gradiente Conjugado x0= soluci´on inicial r0=b−Ax0 s0=r0 for k = 0,1,2,... do αk=rT krk sT kAsk xk+1 =xk+αksk rk+1 =rk−αkAsk βk+1 =rT k+1rk+1 rT krk sk+1 =rk+1 +βk+1sk end for 2.4.2. Gradiente Biconjugado El principal problema que plantea el m´etodo CG es el hecho de que es aplicable ´unicamente a matrices sim´etricas y definidas positivas, haci´endose imposible la ortogonalizaci´on de los vectores residuo de manera eficiente en cualquier otro caso (ver [11] para m´as detalles). Para evitar este problema 18 CAP´ ITULO 2. DIN ´ AMICA DE FLUIDOS COMPUTACIONAL se desarroll´o el algoritmo del gradiente biconjugado (en adelante BiCG por su acr´onimo en ingl´es BiConjugate Gradient). Este m´etodo se basa en la idea de mantener dos subespacios de Krylov cuyos residuos en una determinada iteraci´on kson ortogonales entre s´ı. Concretamente se mantiene el mismo subespacio utilizado en el m´etodo CG Kk(A, r0) = span{r0, Ar0, A2r0, . . . , Ak−1r0} y adem´as se genera otro de la forma Kk(AT,ˆr0) = span{ˆr0, ATˆr0,(AT)2ˆr0,...,(AT)k−1ˆr0} donde ˆr0es un vector aleatorio (generalmente ˆr0=r0) y rT kw= 0 para todo w∈ Kk. Es decir, en una determinada iteraci´on k, cualquier vector del subespacio Kk(y, en concreto, el residuo ˆr0) es ortogonal al residuo rk. Cabe destacar que, puesto que los residuos son ortogonales, las direcciones de b´usqueda generadas ser´an biconjugadas; es decir, ˆsT kAsl= 0 si k6=l. El algoritmo es el siguiente: Algoritmo 2 Gradiente Biconjugado x0= soluci´on inicial r0=b−Ax0 s0=r0 ˆr0=r0 ˆs0= ˆr0 for k = 0,1,2,... do αk=ˆrT krk ˆsT kAsk xk+1 =xk+αksk rk+1 =rk−αkAsk ˆrk+1 = ˆrk−αkATˆsk βk+1 =ˆrT k+1rk+1 ˆrT krk sk+1 =rk+1 +βk+1sk ˆsk+1 = ˆrk+1 +βk+1ˆsk end for Es importante destacar que el uso de BiCG no est´a generalizado debido principalmente al car´acter err´atico del algoritmo en la convergencia. Es por ello que se han desarrollado diferentes optimizaciones del algoritmo, las cuales s´ı que se utilizan de manera generalizada. Concretamente en este trabajo se utiliza el m´etodo BiCGStab(`), aunque su desarrollo e implementaci´on exceden el ´ambito de este trabajo y se pueden encontrar en [12], junto con los principales algoritmos de resoluci´on de sistemas lineales iterativos, tanto estacionarios como no estacionarios. 2.4. SISTEMA DE ECUACIONES LINEALES 19 2.4.3. Residuo m´ınimo generalizado Por ´ultimo, el m´etodo del residuo m´ınimo generalizado (GMRES por su acr´onimo en ingl´es, Generalized Minimum Residue) se basa en la idea de resolver un problema de m´ınimos cuadrados en cada iteraci´on del problema de manera que se minimiza el residuo krkk2=kb−Axkk2, siendo xkla aproximaci´on de la soluci´on exacta x∗en la iteraci´on k. Para conseguir esto, se utiliza el subespacio de Krylov generado por la matriz A y el vector b: Kk(A, b) = span{b, Ab, A2b, . . . , Ak−1b}, y, por tanto, la matriz de Krylov es Kk= [b, Ab, A2b, . . . , Ak−1b]. De esta manera, es posible escribir la aproximaci´on del vector soluci´on como xk=Kkc, donde c∈C. As´ı pues, el residuo a minimizar se convierte en krkk2=kb−Axkk2=kAKkc−bk2. Por otro lado, cabe destacar que las columnas de la matriz de Krylov Kktienden a ser linealmente dependientes conforme aumenta el n´umero de columnas. Es por esto que se utiliza la iteraci´on de Arnoldi con el fin de obtener una base ortonormal Qkpara el subespacio Kkde forma que el vector soluci´on aproximado se reescribe como xk=Qkypara un cierto y∈C, resultando el residuo a minimizar como kb−Axkk2=kAQky−bk2. Teniendo en cuenta la idea principal de Arnoldi AQk=Qk+1 ˆ Hk, donde ˆ Hkes una matriz Hessenberg superior, es posible transformar el residuo a minimizar en: kQk+1 ˆ Hky−bk2=kˆ Hky−QT n+1bk2. Es posible demostrar ([13]) que QT n+1b=kbk2e1, por lo que nuestro problema resultante consiste en minimizar la siguiente norma: kˆ Hky− kbk2e1k2, con xk=Qky. As´ı pues, el m´etodo GMRES consiste en obtener las matrices Qkyˆ Hk mediante la iteraci´on de Arnoldi y, posteriormente, resolver un problema de m´ınimos cuadrados (por ejemplo, mediante la factorizaci´on QR). El algoritmo que implementa este m´etodo es el siguiente: 26 CAP´ ITULO 3. PETSC s´olo el 0) establece todos los valores del vector siguiendo una funci´on lineal x(i) = 10 ∗i. Finalmente, se hace el ensamblado del vector, que consiste en realizar las comunicaciones necesarias para que los valores establecidos por el procesador 0 se distribuyan de manera correcta por todos los procesos. Es importante notar que esta etapa se hace en dos fases, VecAssemblyBegin yVecAssemblyEnd. La utilidad de esto es principalmente agrupar comunicaciones del ensamblado de varios vectores de forma que se utilizara un tiempo de latencia ´unico para todas las comunicaciones, en lugar de uno por cada vector. Esto se consigue agrupando primero todas las llamadas a VecAssemblyBegin seguidas de todas las llamadas a VecAssemblyEnd. Finalmente, y tras realizar todos los c´alculos que sean necesarios, se libera la memoria ocupada por el vector haciendo uso de la funci´on VecDestroy. 3.2. Matrices PETSc proporciona distintas implementaciones para las matrices ya que no existe una ´unica implementaci´on adecuada para todos los tipos de problemas. De esta manera, la librer´ıa proporciona formatos de almacenamiento denso o disperso en sus versiones paralela y secuencial, as´ı como varios formatos especializados. Es importante destacar que desde PETSc se pone mucho ´enfasis en el hecho de que los objetos se definen como un conjunto de operaciones que se pueden realizar con ellos, y no como una implementaci´on concreta. As´ı, una matriz no es un conjunto de n×melementos, sino un objeto que permite el c´alculo de la multiplicaci´on matriz por matriz, la obtenci´on de la transpuesta, etc. Esto se debe al hecho de que la implementaci´on puede (y debe) cambiar en funci´on del problema y de caracter´ısticas tales como el patr´on de dispersi´on de la matriz, pero el resultado de las operaciones sobre ella debe ser siempre el mismo. Tambi´en es de especial inter´es la forma en que PETSc reserva memoria para las matrices dispersas en formato de fila comprimida (CSR por sus siglas en ingl´es, Compressed Sparse Row), conocida en PETSc como formato disperso AIJ. El funcionamiento en PETSc para este tipo de matrices (AIJ y ciertos formatos derivados de ´este), muy comunes en el ´ambito de los algoritmos num´ericos, es que reserva memoria para un n´umero por defecto de elementos por cada fila de manera que, a medida que se construye la matriz, cuando se alcanza el l´ımite se reserva m´as memoria, se hace una copia del contenido actual de la fila y se libera la anterior memoria; este proceso se repite mientras sea necesario hasta la construcci´on completa de la matriz. Esto presenta el problema de que este manejo de memoria y las consiguientes copias de elementos suponen un gasto computacional nada despreciable que puede ser eliminado si se sabe de antemano el n´umero (m´aximo) de elementos que contendr´a la fila m´as densa (o cada una de las filas, si se tiene dicha informaci´on). Este proceso se conoce como preasignaci´on. 3.2. MATRICES 27 As´ı, de manera an´aloga a la creaci´on de vectores, se puede hacer uso de las rutinas MatCreate,MatSetSizes yMatSetFromOptions para crear las matrices en funci´on de los valores especificados en la l´ınea de comandos. Sin embargo, cabe destacar que PETSc proporciona multitud de constructores como MatCreateMPIAIJ,MatCreateMPIDense oMatCreateSeqAIJ con el fin de crear la matriz especificando la implementaci´on concreta que se desea utilizar. MatCreate(PETSC_COMM_WORLD, &m); MatSetSizes(m, PETSC_DECIDE, PETSC_DECIDE, 20, 20); MatSetFromOptions(m); En este punto, es cuando se puede hacer uso de la preasignaci´on mediante la funci´on MatMPIAIJSetPreallocation. MatMPIAIJSetPreallocation(m, 4, PETSC_NULL, 0, PETSC_NULL); Cabe destacar que esta funci´on divide sus argumentos de dos maneras distintas. Por un lado, se entiende la porci´on de la matriz correspondiente a un determinado proceso como la yuxtaposici´on de una submatriz de tama˜no m×m(siendo m, el n´umero de filas que maneja el proceso) conocida como porci´on diagonal (ya que la diagonal de la matriz y la submatriz coinciden) y otra submatriz de tama˜no m×Ncon el resto de la porci´on local de la matriz, tal y como aparece en la figura 3.2. P0 P1 P2 P3 m Nm Figura 3.2: Divisi´on de la matriz para la preasignaci´on. Las regiones en naranja se corresponden con la secci´on diagonal de la submatriz local. Las regiones en azul se corresponden con la secci´on no diagonal. As´ı, el segundo y tercer par´ametro se corresponden con la secci´on diagonal, mientras que el cuarto y quinto se corresponden con el resto. Por 28 CAP´ ITULO 3. PETSC otro lado, se puede realizar la preasignaci´on de manera id´entica para todas las filas (par´ametros segundo y cuarto, enteros) o con valores distintos para cada una de las filas (par´ametros tercero y quinto, vectores de enteros). De esta manera, en el ejemplo anterior se hace la preasignaci´on con 4 elementos no nulos en todas las filas de la secci´on diagonal y ning´un elemento no nulo en el resto de la matriz. Una vez se ha creado la matriz, se pueden establecer los valores de la matriz con MatSetValues y posteriormente realizar el ensamblado con MatAssemblyBegin yMatAssemblyEnd de manera an´aloga al caso de los vectores. Adem´as, obviamente tambi´en se proporcionan funciones para operar con la matriz, tales como el c´alculo de la norma (MatNorm), el escalado (MatScale) y, por supuesto, la liberaci´on de memoria (MatDestroy). Se puede encontrar la documentaci´on de las funciones, as´ı como multitud de ejemplos y tutoriales en [5]. 3.3. Mallas estructuradas En muchos problemas que se modelan mediante ecuaciones de derivadas parciales, es com´un representar el dominio mediante mallas estructuradas rectangulares. En estos casos, una vez se ha realizado la distribuci´on de los datos, suele ser necesaria informaci´on sobre elementos que no se encuentran en local (generalmente elementos contiguos almacenados en nodos vecinos, conocidos como ghost nodes), antes de realizar ciertas operaciones locales. Para ello, PETSc proporciona una estructura de datos espec´ıfica: DA (Distributed Array). Los arrays distribuidos se utilizan junto con los vectores y est´an dise˜nados para comunicar la informaci´on necesaria que no est´a almacenada localmente en mallas estructuradas, tal y como se muestra en la figura 3.3. Es por esto que no se deben utilizar para almacenar matrices o en dominios no estructurados. Al igual que con el resto de objetos de PETSc, los arrays distribuidos se deben crear con DACreate yDASetFromOptions, aunque cabe destacar que la librer´ıa proporciona los constructores DACreate1d,DACreate2d y DACreate3d para crear el objeto especificando el n´umero de dimensiones si este se conoce en tiempo de compilaci´on y no se desea cambiar din´amicamente en tiempo de ejecuci´on. Adem´as, cabe destacar que el uso de los arrays distribuidos est´a fuertemente ligado al uso de vectores. Como se ha comentado anteriormente, las porciones locales del vector deben reservar espacio para los nodos vecinos que ser´an utilizados en el c´alculo. De esta manera, la llamada a la funci´on DAGetLocalVector nos permitir´a obtener un vector local con la memoria necesaria para almacenar todos estos datos. De manera similar, DAGetGlobalVector nos servir´a para obtener un vector que contenga espacio para toda la informaci´on del dominio, distribuida de manera adecuada entre todos los procesos. Por ´ultimo, es importante notar 3.4. SOLVERS LINEALES 29 Figura 3.3: Malla estructurada dividida entre cuatro procesos. Las celdas en gris corresponden a la porci´on local de P0. Las celdas en naranja se corresponden con celdas de procesos distintos a P0 y que se necesitan en este ´ultimo. que cuando los datos globales deben ser distribuidos a las copias locales o viceversa, se debe hacer un manejo especial de los nodos vecinos, ya que ´estos se encuentran replicados en varios nodos. Esto se realiza mediante las rutinas DAGlobalToLocal yDALocalToGlobal que, mediante un par´ametro, permite especificar si en la comunicaci´on se deben insertar los valores de los nodos vecinos, reemplazando los antiguos, o por el contrario, se debe a˜nadir el valor al dato existente. 3.4. Solvers lineales El objeto KSP es el n´ucleo de PETSc ya que proporciona un acceso uniforme a todo el conjunto de m´etodos de resoluci´on de sistemas de ecuaciones lineales, paralelos y secuenciales, directos e iterativos. KSP est´a destinado a resolver sistemas no singulares de la forma Ax =b, donde Aes la matriz de coeficientes, bel vector del lado derecho y xes el vector soluci´on. Para ello, permite entre otros el uso de multitud de m´etodos de Krylov en conjunci´on con diversos precondicionadores, lo cual se ha convertido en el proceso m´as com´un en la resoluci´on de sistemas de ecuaciones lineales de manera iterativa. En cuanto al uso de los objetos KSP, en primer lugar ha de crearse el objeto mediante la llamada KSPCreate. Posteriormente, se deben especificar los 30 CAP´ ITULO 3. PETSC operadores del objeto KSP mediante KSPSetOperators, esto es, la matriz de coeficientes y la matriz a partir de la cual se debe construir el precondicionador que, generalmente, suele coincidir. Tras establecer estas matrices, se debe hacer una llamada a KSPSetFromOptions, que establecer´a los par´ametros, tales como el tipo de m´etodo iterativo, precondicionador y m´ultiples par´ametros de ´estos, todo ello en tiempo de ejecuci´on. A partir de este momento ya es posible resolver el sistema de ecuaciones especificando el vector del lado derecho del sistema y el vector donde se almacenar´a la soluci´on mediante la funci´on KSPSolve. El flujo de ejecuci´on m´as simple para el uso de KSP seKSP ser´ıa el siguiente: ... KSPCreate(PETSC_COMM_WORLD, &ksp); KSPSetOperators(ksp, A, A, DIFFERENT_NONZERO_PATTERN); KSPSetFromOptions(ksp); KSPSolve(ksp, b, x); ... donde Aes la matriz de coeficientes, bel vector de t´erminos independientes y xel vector soluci´on. Adem´as, el ´ultimo par´ametro de la llamada a KSPSetOperators especifica si la nueva matriz de coeficientes tiene o no el mismo patr´on de elementos no nulos que en la anterior llamada a esta misma funci´on. En el ejemplo, puesto que se establece la matriz de coeficientes una sola vez, este par´ametro es irrelevante. Tambi´en cabe destacar que, al establecer por separado la matriz de coeficientes y el vector del lado derecho, se optimiza la creaci´on del precondicionador en caso de que se deseen resolver m´as de un sistema de ecuaciones con la misma matriz de coeficientes y el mismo precondicionador. Por otro lado, las opciones que se le pueden especificar al objeto KSP mediante l´ınea de comandos (y que se aplican al ejecutar la instrucci´on KSPSetFromOptions) son muy numerosas. En primer lugar, es posible especificar el m´etodo de Krylov (u otros) que se desea utilizar en la resoluci´on del sistema con la opci´on -ksp type seguida del m´etodo, pudiendo elegir entre los m´etodos CG (cg), CGS (cgs), BiCG (bicg), GMRES (gmres), Richardson (richardson), Chebychev (chebychev) entre otros. Adem´as, para cada uno de estos m´etodos se pueden especificar opciones particulares, tales como el factor de Richardson (KSPRichardsonSetScale o-ksp richardson scale) o el reinicio de GMRES (KSPGMRESSetRestart o-ksp gmres restart). Por otra parte, tambi´en es posible especificar el precondicionador mediante la opci´on -pc type, pudiendo elegir entre Jacobi (jacobi), Block Jacobi (bjacobi), SOR (sor), LU incompleta (ilu), Cholesky incompleta (icc) entre otros. Para cada uno de estos precondicionadores existen par´ametros concretos que se pueden especificar tales como la omega de SOR (PCSORSetOmega o-pc sor omega). 3.5. SOLVERS NO LINEALES 31 3.5. Solvers no lineales El m´odulo SNES de PETSc proporciona un conjunto de estructuras de datos y operaciones para la resoluci´on de problemas no lineales. As´ı, esta librer´ıa se construye sobre los elementos comentados en apartados anteriores (vectores, matrices, solvers lineales,...) con el fin de proporcionar al usuario una manera flexible y f´acil de establecer los par´ametros del objeto y resolver el sistema de su problema particular. El uso del m´odulo SNES es an´alogo al resto de objetos de PETSc, cre´andolos con SNESCreate,SNESSetFromOptions, resolviendo con SNESSolve y liberando memoria con SNESDestroy. Sin embargo, esta librer´ıa se basa en el m´etodo de Newton para manejar la componente no lineal del sistema. Este m´etodo consiste en la aplicaci´on de la serie de Taylor sobre el sistema no lineal F(x) = 0. As´ı, se utilizan los dos primeros t´erminos de la serie para aproximar la funci´on: F(x)≈F(x0) + F0(x0)(x−x0) (3.1) e igualando a 0 se obtiene la estimaci´on de la nueva soluci´on: x1=x0−F(x0) F0(x0)o, en general, xk=xk−1−F(xk−1) F0(xk−1)(3.2) La soluci´on final se obtiene repitiendo este proceso hasta que el cambio en la soluci´on xk−xk−1es aceptado por el criterio de convergencia. As´ı pues, el principal inconveniente del m´etodo es que se necesita calcular la derivada de la funci´on, conoci´endose esta derivada como el Jacobiano de la funci´on inicial, lo cual plantea dos problemas. Por un lado, la propia evaluaci´on del Jacobiano se suele convertir en la parte que mayor tiempo de computaci´on requiere, mientras que por otro lado, muchos sistemas se construyen a partir de ecuaciones impl´ıcitas o que simplemente son tan complejas que la derivaci´on se hace pr´acticamente imposible. De este modo, tal como se explica en [7], el m´etodo de Newton se ha utilizado muy pocas veces para resolver las ecuaciones de Navier-Stokes de manera directa, demostr´andose en estos casos que, a pesar de que el m´etodo converge en muy pocas iteraciones, el tiempo empleado en el c´alculo del Jacobiano lo convierte en una alternativa peor que otros m´etodos iterativos. Es por esto que en este trabajo se implementa la iteraci´on de Picard para el problema no lineal y no se hace uso de la librer´ıa SNES de PETSc. 3.6. Integradores temporales La librer´ıa TS (por sus siglas en ingl´es, Time Steppers) proporciona un entorno de trabajo para la programaci´on de soluciones escalables, tanto en 32 CAP´ ITULO 3. PETSC problemas de ecuaciones de derivadas parciales dependientes del tiempo, como en problemas estacionarios con del tiempo. De manera an´aloga al caso de los solvers no lineales, en el problema particular que se pretende resolver en este trabajo, la integraci´on temporal es compleja y se requiere de una gran flexibilidad, por lo que el c´odigo cuenta con una implementaci´on propia y no se hace uso de la librer´ıa TS. Cap´ıtulo 4 MICSc En un principio, el proyecto se inicia con la idea de desarrollar por completo una nueva aplicaci´on para la resoluci´on num´erica de las ecuaciones de Navier-Stokes. Sin embargo, poco tiempo despu´es, y tras revisar las distintas aplicaciones de la librer´ıa PETSc, se fija la atenci´on en una tesis realizada en la Universidad de Zaragoza por A. Cubero bajo la tutela de N. Fueyo [14]. Este c´odigo se centra exactamente en el objetivo principal del proyecto de investigaci´on RHELES y finalmente, tras contactar con los autores, se decide aprovechar el trabajo ya realizado y se adopta MICSc como punto de partida para el desarrollo del proyecto. As´ı pues, MICSc es el nombre del c´odigo escrito en C, de aproximadamente 15000 l´ıneas en ficheros fuente (incluyendo documentaci´on y espacios), que hace uso de PETSc para la resoluci´on acoplada e impl´ıcita de las ecuaciones de Navier-Stokes. En este cap´ıtulo se detallan las principales caracter´ısticas de este c´odigo, muchas de las cuales se han mantenido mayormente inalteradas y son sobre las que se han aplicado correcciones y nuevas funcionalidades (comentadas en el siguiente apartado) con el fin de alcanzar los objetivos del proyecto. De esta manera, en primer lugar trataremos la forma en que se abordan los elementos m´as problem´aticos para la implementaci´on y las decisiones que se han tomado y, al final del cap´ıtulo, entraremos con m´as detalle en la implementaci´on concreta de la soluci´on. 4.1. Discretizaci´on La discretizaci´on en MICSc se realiza siguiendo el modelo de vol´umenes finitos. Este m´etodo se fundamenta en la divisi´on estructurada del dominio en un conjunto de celdas sobre las que se aplican individualmente las ecuaciones de gobierno del flujo. Sin embargo, existen dos aproximaciones distintas a la hora de ubicar los valores del problema sobre la malla: mallas decaladas y mallas colocalizadas (se puede ver una comparativa de ambas en [15]). La principal diferencia es que en las mallas colocalizadas el valor 33 34 CAP´ ITULO 4. MICSC x2, y xδ vn MI u e MI s p EP N S W e n s w x∆ x1, x vu Figura 4.1: Notaci´on de las celdas, en may´uscula los centros de las celdas y en min´uscula las caras (izquierda). Distribuci´on de las variables, sub´ındice indicando el vecino, super´ındice MI indicando la Interpolaci´on del Momento de una componente de la velocidad. de todas las inc´ognitas se almacena en el centro de la celda, lo cual puede llevar al desacoplamiento de la presi´on (de la que s´olo aparece el gradiente) y la velocidad y, por tanto, a resultados con poco sentido f´ısico. Por contra, en las mallas decaladas se transportan las componentes de la velocidad a las caras de las celdas, manteniendo los valores escalares en el centro. En el caso de MICSc, se hace uso de la Interpolaci´on del Momento [16] (MI por sus siglas en ingl´es, Moment Interpolation), que consiste en almacenar tanto las componentes de la velocidad como los escalares en las celdas, y definiendo unas nuevas variables sobre las caras (velocidades de convecci´on) que se obtienen mediante una interpolaci´on concreta de las componentes de la velocidad. La Figura 4.1 ilustra la disposici´on de estas variables y su notaci´on. Por otro lado, una vez decidida la implementaci´on de la malla y la distribuci´on de las variables sobre ella, discretizaremos las ecuaciones. Es posible demostrar que cualquier ecuaci´on de conservaci´on sobre una variable φes posible expresarla de manera semi-discretizada como la suma de los t´erminos temporal, convectivo, difusivo y fuente de la siguiente manera: ∂ρφ ∂t PVP+X nb(P) mnbφnb −X nb(P) Γφ nbAnb ∂φ ∂xnb nb =Sφ PVP donde VPes el volumen de la celda P; los sumatorios sobre nb(P) recorren todas las caras de la celda; mnb es el flujo a trav´es de la cara nb, definido como mnb =ρnbAnbvnb donde vnb es la componente de la velocidad normal 4.2. ENTRADA/SALIDA 35 a la cara nb yAnb =Aˆn·ˆe; Γφ nb es un coeficiente de difusi´on para la variable φ;xnb representa la direcci´on normal a la cara nb ySφ Pengloba los t´erminos fuente (o sumidero) de φ. Llegados a este punto, la discusi´on sobre las diferentes maneras de discretizar los distintos t´erminos se puede encontrar en [14] y excede los l´ımites de este trabajo. ´ Unicamente destacar que para el t´ermino temporal se utilizan esquemas temporales impl´ıcitos (Euler, Adams-Moulton, ...) y para el t´ermino convectivo se utilizan esquemas convectivos de alto orden (SMART, QUICK, ...). Sin embargo, s´ı que es importante destacar el hecho de que, una vez discretizados estos t´erminos, las ecuaciones de conservaci´on del momento son realmente apropiadas para despejar las componentes de la velocidad. No obstante, la ecuaci´on restante (conservaci´on de la masa) no resulta adecuada para la presi´on, puesto que no aparece dicha variable y el sistema algebraico resultante presenta un 0 en la diagonal. Nuevamente, existen diversas aproximaciones y multitud de estudios sobre la soluci´on de este problema que exceden el ´ambito del proyecto. Para el caso concreto de MICSc, la soluci´on pasa por derivar una ecuaci´on de tipo Poisson a partir de la ecuaci´on de continuidad. Para m´as detalles matem´aticos y de implementaci´on sobre esta soluci´on, ver [14] [17]. 4.2. Entrada/Salida Otro aspecto importante de la aplicaci´on es la forma en que recibe los datos y produce los resultados. Al margen de las opciones que proporciona PETSc tanto para la entrada de par´ametros como para la salida de resultados, es necesario contar con otro tipo de datos y formatos de salida que se especifican en los siguientes apartados. 4.2.1. Entrada Uno de los ficheros de entrada m´as importantes es el que contiene la malla que define los vol´umenes finitos que se utilizar´an en la resoluci´on del problema. Aunque actualmente se sigue un formato propio para este fichero, se consideran otras alternativas como trabajo futuro (ver 7.2). En cuanto al formato propio de este fichero, actualmente se define de la siguiente manera: <{0,1}> <num_divisiones_eje_X> <num_divisiones_eje_Y> <num_divisiones_eje_Z> <lista_coords_eje_X> <lista_coords_eje_Y> 42 CAP´ ITULO 4. MICSC todo el desarrollo debido a que almacena los par´ametros de entrada del problema y ´estos se han ido modificando, tanto por el hecho de a˜nadir nueva funcionalidad a la aplicaci´on como por la extracci´on de m´ultiples par´ametros que en el inicio se especificaban dentro del c´odigo fuente y se extrajeron al fichero de par´ametros de entrada. St InputBcond: typedef struct St_InputBcond { PetscInt x, y, z, nx, ny, nz, ineighb; BcondType type; PetscReal porosity; struct St_InputBcond *next; } St_InputBcond; Contiene la definici´on de una condici´on de contorno tal y como se especifica en el fichero de entrada. Esta estructura ´unicamente tiene utilidad para almacenar los valores temporalmente antes de crear las estructuras St Bcond asociadas. Sus atributos son el punto inicial (x, y,z) de la regi´on sobre la que se define la condici´on de contorno, as´ı como el n´umero de celdas de dicha regi´on en cada direcci´on (nx, ny,nz). Tambi´en contiene el tipo de condici´on (type), la porosidad (porosity) y un puntero para crear listas (next). St Interp: typedef struct { PetscReal **linear, **upwind; PetscInt **icell_upw, **icell_uupw, **icell_dow; } St_Interp; Esta estructura de datos es ´util para el c´alculo de los coeficientes de interpolaci´on del esquema convectivo. Para cada dimensi´on y cada celda contiene datos como el ´ındice local de las celdas upwind (icell upw), up-upwind (icell uupw) y downstream (icell dow). Tambi´en incluye los coeficientes para la interpolaci´on tanto lineal (linear) como upwind (upwind). St Local: typedef struct { PetscMPIInt rank; MeshRegion meshRegion, ghostMeshRegion; PetscInt **icell_neighb; PetscReal **porosity; St_Porosity *FirstPorosity, *LastPorosity; } St_Local; 4.4. ESTRUCTURACI ´ ON DEL C ´ ODIGO 43 Par´ametros locales del problema. Contiene el ´ındice del procesador (rank) e informaci´on sobre la porci´on de malla que maneja, tanto considerando los nodos ghost (ghostMeshRegion), como sin considerarlos (meshRegion). Tambi´en se almacena una lista de porosidades (FirstPorosity, LastPorosity) que se aplican sobre la porci´on local de la malla. Por ´ultimo, para cada celda, y en funci´on del n´umero de dimensiones del problema, se almacenan los ´ındices de sus vecinos (icell neighb) y el valor num´erico de su porosidad (porosity). St MeshRegion: typedef struct { PetscInt nCells[3]; PetscInt lowerCorner[3]; PetscInt upperCorner[3]; PetscInt totalCells; } St_MeshRegion; Define una regi´on rectangular tridimensional de la malla, especificando el punto origen (lowerCorner), el punto final (upperCorner), as´ı como el n´umero de celdas que contiene dicha regi´on en cada direcci´on (nCells) y en total (totalCells). St Patch: typedef struct St_Patch { PetscInt id; PetscInt x, y, z, nx, ny, nz; DA da; MPI_Comm comm; struct St_Patch *Next; } St_Patch; Esta estructura de datos contiene un identificador (id), as´ı como las dimensiones globales de una determinada regi´on de la malla (x,y,z, nx,ny,nz). Adem´as, contiene un objeto DA y el comunicador asociado (comm) para manejar las propiedades que se definan sobre esta regi´on de la malla. Por ´ultimo, proporciona un puntero para construir listas (Next). St Porosity: typedef struct St_Porosity { PetscInt id, ineighb; St_Patch *Patch; PetscReal Cons; struct St_Porosity *Next; } St_Porosity; 44 CAP´ ITULO 4. MICSC Esta estructura almacena par´ametros ´utiles para relajar o limitar flujos a trav´es de las caras (con un factor entre 0 y 1). Concretamente, contiene un identificador (id), el ´ındice de la cara sobre la que se aplica el factor (ineighb), el conjunto de celdas sobre la que se aplica la porosidad (Patch) y el factor por el que se multiplica el flujo (Cons). Cabe destacar tambi´en el hecho de que cada vez que se crea una nueva porosidad, ´esta se a˜nade a una lista de porosidades locales que se puede recorrer haciendo uso del atributo Next. St Prop: typedef struct St_Prop { PetscInt id; St_Patch *Patch; ktypeProp ktype; ktypeIp kip; PetscErrorCode (*FluxLimiter)(PetscInt icell, PetscInt ineighb, PetscReal upwdelta, PetscReal dowdelta, PetscReal *value); PetscReal linrlx; PetscReal Cons; PetscErrorCode (*Funct)(struct St_Prop *Prop); PetscReal *array; Vec instantMeans, instantSums, varianceSums, values; PetscInt numReferences; struct St_Prop *Next; } St_Prop; Propiedad (constante o variable) asociada a una regi´on de la malla (o patch). Entre sus atributos encontramos un identificador (id), la regi´on de la malla sobre la que se aplica la propiedad (Patch), y el tipo de propiedad (ktype), que define si la propiedad es constante (CONS) o depende de la geometr´ıa (GEOM), de las inc´ognitas (VARI) o de los t´erminos transitorios (TRANS). Adem´as, contiene el tipo de interpolaci´on (kip), que define la funci´on utilizada para el limitador de flujo (FluxLimiter) y un valor de relajaci´on (linrlx). Por ´ultimo, si el tipo de propiedad es constante, el atributo Cons almacenar´a este dato constante y, en caso contrario har´a uso de la funci´on Funct para calcular los datos. Tambi´en contiene objetos Vec donde se almacenar´an los valores instant´aneos de la propiedad (values), as´ı como valores para el c´alculo de medias (instantMeans,instantSums yvarianceSums, ver 5.2.2). Es posible acceder al almacenamiento del vector de valores instant´aneos mediante el atributo array. Puesto que una propiedad puede estar asociada a distintos objetos de la aplicaci´on (principalmente a la gamma de las variables), existe un contador con el que se 4.4. ESTRUCTURACI ´ ON DEL C ´ ODIGO 45 determina el n´umero de referencias a esta propiedad(numReferences). Adem´as, cada vez que se crea una propiedad, se a˜nade a una lista de propiedades locales, iterables mediante Next. St Sys: typedef struct { Vec Phi, B; Mat A; DA da; KSP solver; PetscInt *ltog; PetscInt ndof; PetscInt ivar; PetscReal resnorm, resref, corrnorm, corref; } St_Sys; Estructura para la creaci´on y resoluci´on del sistema de ecuaciones distribuido. Contiene el vector inc´ognita (Phi), el vector del lado derecho (B), la matriz de coeficientes (A), el array distribuido DA de PETSc (da), el solver lineal (solver). Adem´as, contiene un vector de enteros para emparejar los ´ındices de las celdas del contexto local al contexto global (ltog). Tambi´en almacenan los grados de libertad del sistema (ndof), que toma el valor 1 para el sistema segregado y el n´umero de variables acopladas en el sistema acoplado. As´ı, sabiendo el grado de libertad y el ´ındice de la primera variable que resuelve el sistema (ivar) se puede conocer el conjunto de variables que resuelve el sistema. Por ´ultimo, tambi´en contiene valores de referencia (resref, corref) para la normalizaci´on de las tolerancias relativa y absoluta (resnorm,corrnorm). St Trans: typedef struct { TransientOrder norder; PetscReal time; PetscInt idt; } St_Trans; Par´ametros espec´ıficos de un flujo no estacionario, como el orden del esquema temporal (norder), que puede ser Euler de primer orden (EUL1), Euler de segundo orden (EUL2), Adams-Moulton de segundo orden (ADM2) o Adams-Moulton de tercer orden (ADM3). Adem´as, tambi´en contiene el instante de tiempo actual a lo largo de toda la ejecuci´on (time) y el n´umero de pasos de tiempo efectuados (idt). St Var: typedef struct St_Var { 46 CAP´ ITULO 4. MICSC char name[128]; PetscInt index; ktypeIp kip; PetscErrorCode (*FluxLimiter)(PetscInt icell, PetscInt ineighb, PetscReal upwdelta, PetscReal dowdelta, PetscReal *value); PetscReal valmin, valmax, linrlx, fdtrlx, resref, corref, resnorm, corrnorm; St_Prop *Valini, *Gamma; Vec locPhi; PetscReal *Array; struct St_Var *next; } St_Var; Contiene par´ametros de una variable dependiente (o inc´ognita). Entre los datos que incluye podemos encontrar el nombre (name), el ´ındice (index), as´ı como el vector de PETSc con los valores locales de la variable (locPhi) y un puntero a los valores almacenados en dicho objeto (Array). Adem´as, tambi´en contiene el tipo de interpolaci´on (kip) que, de manera an´aloga al caso de St Prop, define el limitador de flujo (FluxLimiter). Junto con todo esto, incluye los valores m´ınimo y m´aximo para la variable (valmin,valmax), as´ı como coeficientes de relajaci´on (linrlx,fdtrlx), valores de referencia (resref,corref) para la normalizaci´on de las tolerancias relativa y absoluta (resnorm, corrnorm) y dos propiedades utilizadas para el c´alculo de los valores iniciales (Valini) y el coeficiente de difusi´on (Gamma). Por ´ultimo, cada vez que se crea una variable nueva antes del inicio de la resoluci´on del problema, ´esta se a˜nade a una lista enlazada haciendo uso del atributo next. 4.4.2. Instancias y relaci´on de las estructuras de datos A pesar de que a partir de las definiciones de las estructuras de datos se deja intuir el n´umero de instancias que se crear´an de cada una de ellas, en este apartado especificaremos de manera concreta dichas instancias: St Bcond: Las condiciones de contorno se especifican en el fichero de entrada y cada una de las condiciones de contorno se convierte en una instancia de la estructura St Bcond. St BcondApplication: Como se ha comentado anteriormente, cada condici´on de contorno contiene diversas aplicaciones (applicationList) que introducen valores en el sistema de ecuaciones para cada ecuaci´on y variable particular. As´ı pues existir´a una instancia distinta para cada una de estas contribuciones y este n´umero de instancias depender´a de 4.4. ESTRUCTURACI ´ ON DEL C ´ ODIGO 47 la condici´on de contorno concreta, por lo que se deber´ıa especificar para cada caso particular. St Geom: Los datos que contiene esta estructura de datos son propios de cada proceso, por lo que existir´a una instancia por cada uno de ellos que contendr´a la informaci´on de la porci´on local de la malla. St Global: Puesto que la informaci´on almacenada en esta estructura de datos es global, existir´a un instancia en cada proceso que contendr´a dicha informaci´on replicada en cada uno de ellos. St Grid: An´alogamente al caso anterior, los datos que se incluyen en esta estructura de datos son globales, por lo que habr´a una instancia id´entica en cada uno de los procesos. St Input: Nuevamente, la informaci´on de St Input es global y habr´a una misma instancia replicada en todos y cada uno de los diferentes procesos. St InputBcond: Puesto que esta estructura ´unicamente se utiliza para leer los datos de entrada del fichero para posteriormente crear las estructuras St Bcond, existir´a una instancia id´entica en todos los procesos por cada condici´on de contorno especificada por el usuario en el fichero de entrada. St Interp: En este caso, la estructura de datos contiene varios vectores cuyo tama˜no depende del n´umero de celdas en cada nodo local. Por tanto, existir´a una instancia en cada uno de los procesos, pero la informaci´on que almacene cada instancia ser´a diferente, puesto que las celdas que se manejan en cada nodo son distintas. St Local: Este caso es muy similar al de St Interp. Existir´a una instancia en cada proceso y dicha instancia incluir´a la informaci´on sobre la regi´on local de la malla que se maneja, as´ı como listas de propiedades, condiciones de contorno y porosidades que se aplican en esa regi´on y dem´as informaci´on local que, puesta en com´un con la informaci´on del resto de procesos, debe reunir todos los datos del dominio completo. St MeshRegion: No existen instancias independientes de esta estructura de datos, sino que se encuentran siempre ligadas como atributos de las estructuras St Local ySt Patch. St Patch: Esta estructura de datos contiene datos de la regi´on de la malla (patch). Es por esto que, por cada regi´on que se deba definir, existir´a una instancia en cada proceso con datos locales e informaci´on global replicada en todos y cada uno de ellos. Cabe destacar tambi´en que, salvo la regi´on que define el domino completo, todas las instancias 48 CAP´ ITULO 4. MICSC St_Global comm: MPI_Comm NPatch, NBcond, NPorosity, NProp: PetsInt FirstBcond, LastBcond: *St_Bcond FirstPatch, LastPatch: *St_Patch FirstProp, LastProp: *St_Prop St_Bcond id, ineighb: PetscInt Patch: *St_Patch bCondType: BcondType applicationList: *St_BcondApplication Next: *St_Bcond St_BcondApplicaton ivar, ieq: PetscInt CoeffA, CoeffB: *St_Prop next: *St_BcondApplication wallDistance: *PetscReal nearestWall: *St_Bcond St_Patch id: PetscInt local, ghost, global: MeshRegion Next: *St_Patch St_Input log_level: PetscInt ... firstBcond, lastBcond: *St_InputBcond St_InputBcond x, y, z, nx, ny, nz: PetscInt type: BcondType porosity: PetscReal next: *St_InputBcond St_Geom Surface: *PetscReal[3] Edge: *PetscReal[3] Volume: *PetscReal wallDistance: *PetscReal nearestWall: *St_Bcond St_Trans norder: TransientOrder time: PetscReal idt: PetscInt St_Prop id: PetscInt Patch: *St_Patch ktype: ktypeProp kip: ktypeIp FluxLimiter: <function> ndof: PetscInt linrlx, Cons: PetscReal Funct: <function> Array: **PetscReal Next: *St_Prop St_Grid Nnode: PetscInt[3] Xnode: *PetscReal[3] ktype: ktypeGrid St_Local rank: PetscMPIInt meshRegion, ghostMeshRegion: St_MeshRegion icell_neighb: **PetscInt porosity: **PetscReal FirstPorosity, LastPorosity: *St_Porosity St_Porosity id, ineighb: PetscInt Patch: *St_Patch Cons: PetscReal Next: *St_Porosity St_MeshRegion nCells: PetscInt[3] lowerCorner: PetscInt[3] upperCorner: PetscInt[3] totalCells: PetscInt St_Interp linear, upwind: **PetscReal icell_upw, icell_uupw, icell_dow: **PetscInt St_Var name: char[128] index: PetscInt ktypeIp: kip FluxLimiter: <function> valmin, valmax, linrlx, fdtrlx: PetscReal resref, corref, resnorm, corrnorm: PetscReal Valini, Gamma: *St_Prop locPhi: Vec Array: *PetscReal next: *St_Var St_Sys Phi, B: Vec A: Mat da: DA solver: KSP ltog: *PetscInt ndof, ivar: PetscInt resnorm, resref, corrnorm, corref: PetscReal Figura 4.2: Relaciones entre las estructuras de datos en MICSc. Por simplicidad se obvian las relaciones recursivas. 4.4. ESTRUCTURACI ´ ON DEL C ´ ODIGO 49 de St Patch est´an asociadas a otras instancias de estructuras de datos diferentes y, por tanto, el n´umero total de instancias que se crear´an de esta estructura depender´a del n´umero de variables, condiciones de contorno y porosidades que se definan para cada problema particular. St Porosity: Las porosidades son estructuras de datos que se utilizan para modificar el flujo a trav´es de las caras de las celdas y est´an siempre asociadas a condiciones de contorno. As´ı pues, existir´a una porosidad por cada condici´on de contorno especificada en el fichero de entrada. Sin embargo, puesto que las condiciones de contorno y, por tanto, las porosidades se aplican sobre una regi´on del dominio concreta, ´unicamente se crear´an instancias en los procesos que contengan la regi´on de aplicaci´on (completa o parcialmente). St Prop: Esta estructura de datos se utiliza para evaluar una cierta propiedad en una regi´on del dominio y almacenar los resultados obtenidos para cada celda local en un vector. Por este motivo, y de manera an´aloga al caso anterior, s´olo se crear´an instancias de St Prop en los procesos que manejen celdas donde se aplique la propiedad de manera completa o parcial, almacenando el vector de manera distribuida en el segundo caso. Por otra parte, existen ciertas propiedades que se eval´uan y se introducen en el sistema de ecuaciones de manera independiente, tales como los t´erminos convectivo (Conv) y difusivo (Diff) o la densidad (Rho). Sin embargo, la mayor parte de las propiedades se encuentran asociadas a condiciones de contorno (CoeffA,CoeffB) o a variables (Valini,Gamma) y el n´umero de instancias que se crear´an depender´a del n´umero de condiciones de contorno y variables que se creen y, por tanto, de cada problema particular. St Sys: Esta estructura de datos contiene principalmente la informaci´on del sistema de ecuaciones: matriz de coeficientes, vector del lado derecho, vector de inc´ognitas y el contexto del array distribuido (DA). Todos ellos son objetos de PETSc que se manejar´an de manera distribuida en todos los procesos haciendo uso de la librer´ıa de manera transparente. Adem´as, tambi´en contiene ciertos par´ametros globales de problema, tales como valores num´ericos globales del problema tales como valores de referencia (resref,corref) para la normalizaci´on de las tolerancias relativa y absoluta (resnorm,corrnorm) que se replicar´an en todos los nodos. Por ´ultimo, es destacable tambi´en que existir´an dos instancias de esta estructura de datos en cada proceso: una de ellas se encargar´a de la resoluci´on del sistema de ecuaciones acoplado y la otra, del segregado. St Trans: En este caso, los datos que incluye la estructura de datos son simplemente valores num´ericos para el manejo de flujos no estacionarios 50 CAP´ ITULO 4. MICSC que se replicar´an en todos los procesos, creando una ´unica instancia id´entica en todos y cada uno de los procesos. St Var: Esta estructura de datos contiene varios par´ametros (valores m´ınimo y m´aximo, relajaciones, ...) que se replicar´an en todos los procesos, junto con un vector que contendr´a el vector soluci´on para las celdas que se manejen en cada uno de los procesos, as´ı como las propiedades Gamma yValini que se distribuir´an tal y como se ha detallado anteriormente para la estructura St Prop. Adem´as, existir´a una instancia con dichos datos por cada una de las variables que se defina en el problema (velocidades, presi´on, temperatura, ...) Conociendo ya la definici´on de las estructuras y el n´umero de instancias que se generan de cada una de ellas, pasamos a resumir y esquematizar las relaciones que existen entre cada una de estas estructuras de datos. En un principio, es importante determinar las estructuras de datos que ´unicamente contienen valores num´ericos, objetos de PETSc y dem´as datos que no son estructuras de datos propias de MICSc. ´ Estas son: St Grid,St Interp, St Sys ySt Trans. Adem´as, cabe destacar el hecho de que tanto St Global como St Local est´an relacionados con las propiedades, condiciones de contorno y/o porosidades, pero ´unicamente con un puntero al inicio de cada una de las listas de objetos. Por otro lado, tanto St Bcond como St Porosity ySt Prop tienen una referencia a St Patch, que define la regi´on de la malla (patch) donde se aplican las condiciones de contorno, porosidades y propiedades respectivamente. Por ´ultimo, tanto St BcondApplication como St Var tienen referencias aSt Prop. Es importante notar que las condiciones de contorno y las variables tienen relaci´on con las regiones de la malla y con propiedades que, a su vez, tambi´en se relacionan con estas regiones. En este caso, debido a la transitividad, se ha cuidado que la implementaci´on de las diferentes referencias a regiones de la malla apunten exactamente a la misma instancia, por lo que no existe ninguna instancia que se replique de manera innecesaria en ning´un nodo. En la figura 4.2 se muestran tanto las estructuras de datos que no dependen de ninguna otra estructura de datos definida en MICSc como aquellas que si que tienen dependencias entre s´ı. 4.4.3. Flujo de ejecuci´on Como se ha comentado anteriormente, MICSc ´unicamente hace uso de los solvers lineales de PETSc e implementa su propia versi´on del m´etodo de Picard para las iteraciones exteriores y la integraci´on temporal, por lo que esto se ver´a reflejado en el flujo de ejecuci´on de la simulaci´on. 4.4. ESTRUCTURACI ´ ON DEL C ´ ODIGO 51 B´asicamente, este flujo consiste en una primera secci´on de lectura de par´ametros de entrada, creaci´on de estructuras e inicializaci´on de valores. Posteriormente, en la parte central se encuentra la resoluci´on del sistema de ecuaciones. Esta secci´on contiene, en primer lugar, el bucle correspondiente a las iteraciones temporales, que al inicio del bucle se encarga de evaluar las propiedades dependientes del tiempo en cada iteraci´on y, al final, actualizar el siguiente paso de tiempo en caso de que la simulaci´on lo requiera, o haciendo falsa la condici´on del bucle en caso contrario. Dentro de este bucle se encuentra anidado otro bucle que se corresponde con las iteraciones del m´etodo de Picard. Aqu´ı, en cada iteraci´on se eval´uan las propiedades que dependen de las variables y se resuelve el sistema acoplado (en caso de que haya variables acopladas) y tantos sistemas segregados como variables no acopladas existan. La resoluci´on de este sistema consiste en la construcci´on de la matriz de coeficientes y el vector del lado derecho, introduciendo incrementalmente las contribuciones de las condiciones de contorno, la gravedad (en caso de que sea relevante), etc. Finalmente, dentro del bucle del m´etodo de Picard, y tras la resoluci´on de todos los sistemas, se limitan los valores de las variables en caso de que hayan sobrepasado los l´ımites m´ınimo y/o m´aximo. Y, por ´ultimo, la finalizaci´on de los objetos de PETSc, liberaci´on de memoria y dem´as tareas de terminaci´on del programa. Tanto en el Algoritmo 4 como en la Figura 4.3 se muestra el flujo de ejecuci´on. 58 CAP´ ITULO 5. DESARROLLOS MICScVarSetInitValues(PetscInt index, ktypeProp initValuesType, ktypeIp initValuesInterpType, PetscReal initValuesConst, PetscErrorCode (initValuesFunct*) (St_Prop *prop)): Se define de manera an´aloga a MICScVarSetGamma para los valores iniciales de la variable en cada celda. MICScVarSetInterpolation(PetscInt index, ktypeIp interpType): Modifica el tipo de interpolaci´on de la variable. MICScVarSetLimits(PetscInt index, PetscReal minValue, PetscReal maxValue): Modifica los valores m´ınimo y m´aximo de la variable. MICScVarSetReference(PetscInt index, PetscReal resRef, PetscReal corrRef): Modifica los valores de referencia de la variable. MICScVarSetRelaxation(PetscInt index, PetscReal linRelax, PetscReal fdtRelax): Modifica los valores de relajaci´on de la variable. MICScSetInletFunction(PetscErrorCode (*function)(PetscInt*, PetscReal*)): Establece la funci´on que se debe utilizar para calcular el flujo de entrada, en caso de que se defina tal condici´on de contorno. As´ı, un posible flujo b´asico de la ejecuci´on de un caso con MICSc queda de la siguiente manera: MICScInitialize(); MICScVarSetLimits(U, 20, 20); ... MICScSolve(); MICScFinalize(); 5.1.2. Tratamiento din´amico de las variables Como consecuencia de la refactorizaci´on comentada en el apartado anterior, se plante´o la posibilidad de que el usuario pudiera definir nuevas variables. Sin embargo, debido a que esto no se hizo necesario para los casos 5.1. MANTENIMIENTO Y OPTIMIZACI ´ ON DEL C ´ ODIGO 59 de prueba que se utilizaron en el proyecto, no se a˜nadieron estas funciones al conjunto de funciones disponibles por el usuario y consta como trabajo futuro. No obstante, s´ı que es posible ejecutar un cierto caso contemplando la temperatura como variable o no (-use t). Por otra parte, existen ciertas secciones del c´odigo que necesitan conocer el n´umero total de variables con el fin de reservar memoria en funci´on del n´umero de estas variables, tales como diversas condiciones de contorno, propiedades asociadas a las derivadas temporales de las variables o la misma creaci´on del vector que contiene a las variables. En definitiva, es necesario conocer el n´umero total de variables en el momento de la inicializaci´on. Una primera soluci´on que se contempl´o pasaba por incluir un par´ametro de entrada en el programa (-nvar) para que el usuario determinara en el mismo inicio de la ejecuci´on el n´umero final de variables que se iban a necesitar. El problema que planteaba esta soluci´on es que daba lugar a inconsistencias si el usuario determinaba un n´umero de variables por l´ınea de comandos y, posteriormente, el n´umero real era distinto. De esta manera, tal y como se ha dejado intuir en 4.4.1, la soluci´on final consiste en crear las variables como una lista enlazada mediante un puntero a la siguiente variable (next). La siguiente secci´on de c´odigo muestra la creaci´on de variables con esta nueva implementaci´on: St_Var pointer, newVariable; PetscNew(St_Var, &newVariable); if (firstVar == PETSC_NULL) { firstVar = newVariable; } else { pointer = firstVar; while (pointer->next != PETSC_NULL) { pointer = pointer->next; } pointer->next = newVariable; } newVariable->next = PETSC_NULL; newVariable->index = num_vars; ... <establecer atributos> ... num_vars++; donde firstVar es una variable global de tipo St Var que apunta al primer elemento de la lista. 60 CAP´ ITULO 5. DESARROLLOS P2 P3 P1P0 Figura 5.1: Distribuci´on de un patch. En color azul el dominio completo, en color salm´on el patch, dividido entre los procesos 0 y 1. 5.1.3. Refactorizaci´on de los patches En la versi´on inicial de MICSc que se recogi´o al principio del proyecto, los patches se implementaban mediante tres instancias de la estructura St MeshRegion: para la regi´on global del patch y la subregi´on local con y sin nodos ghost. Sin embargo, se comprob´o que esta implementaci´on proporcionaba exactamente la misma funcionalidad que los objetos DA de PETSc. Adem´as, tambi´en se daba el hecho de que en el inicio las propiedades se evaluaban de manera local (incluyendo nodos vecinos) y aislada en cada proceso, sin producirse comunicaci´on entre procesos, lo cual pod´ıa ser la causa de cierto comportamiento paralelo de divergencia para instantes de tiempo avanzados en la simulaci´on. Por ´ultimo, y con el fin de simplificar la implementaci´on de muchas secciones del c´odigo, se deseaba que la evaluaci´on de las propiedades fuera exclusivamente local (sin nodos vecinos), y que se produjera comunicaci´on entre procesos despu´es de evaluar cada una de las propiedades. Por todo esto, se decidi´o sustituir la implementaci´on inicial por una implementaci´on con un objeto DA de PETSc. Esta implementaci´on planteaba numerosas dificultades. La primera de ellas es que el array distribuido se deb´ıa organizar de la misma manera que se organizaba el resto del dominio; es decir, dado un cierto patch, ´este deb´ıa utilizar las mismas divisiones que las usadas en el dominio completo, pudiendo ser posible que la regi´on se dividiera ´unicamente entre un subconjunto de procesos (ver Figura 5.1). De esta forma, lo primero que se debe hacer es calcular la porci´on local que corresponde a cada proceso: 5.1. MANTENIMIENTO Y OPTIMIZACI ´ ON DEL C ´ ODIGO 61 *r = Local.meshRegion; IsPatchInProc(x, nx, y, ny, z, nz, &flag); if (flag) { overlap.lowerCorner[XDIR] = MAX(x, r->lowerCorner[XDIR]); overlap.lowerCorner[YDIR] = MAX(y, r->lowerCorner[YDIR]); overlap.lowerCorner[ZDIR] = MAX(z, r->lowerCorner[ZDIR]); overlap.upperCorner[XDIR] = MIN(x + nx, r->lowerCorner[XDIR] + r->nCells[XDIR]) - 1; overlap.upperCorner[YDIR] = MIN(y + ny, r->lowerCorner[YDIR] + r->nCells[YDIR]) - 1; overlap.upperCorner[ZDIR] = MIN(z + nz, r->lowerCorner[ZDIR] + r->nCells[ZDIR]) - 1; overlap.nCells[XDIR] = overlap.upperCorner[XDIR] - overlap.lowerCorner[XDIR] + 1; overlap.nCells[YDIR] = overlap.upperCorner[YDIR] - overlap.lowerCorner[YDIR] + 1; overlap.nCells[ZDIR] = overlap.upperCorner[ZDIR] - overlap.lowerCorner[ZDIR] + 1; overlap.totalCells = overlap.nCells[XDIR] * overlap.nCells[YDIR] * overlap.nCells[ZDIR]; } else { overlap.totalCells = 0; } donde x,y,z,nx,ny,nz representan el patch. As´ı, se determina si ´este se debe distribuir en el proceso local (IsPatchInProc), marc´andose el n´umero total de celdas como 0 en caso negativo. En caso de que el patch se solape con el subdominio local, se asignan los atributos del objeto overlap, de tipo St MeshRegion, que contienen la intersecci´on del subdominio local con el patch. Posteriormente, se deben distribuir estas intersecciones a todos los procesos: MPI_Type_struct(4, lengths, displs, types, &meshRegionType); MPI_Type_commit(&meshRegionType); MPI_Allgather(&overlap, 1, meshRegionType, overlaps, 1, meshRegionType, PETSC_COMM_WORLD); donde lengths,displs ytypes son, respectivamente, las longitudes, el desplazamiento en bytes y los tipos de los atributos de la estructura de datos St MeshRegion. De esta manera, se crea un nuevo tipo de datos MPI que se corresponde con dicha estructura (meshRegionType) y se distribuyen los objetos calculados anteriormente a todos los procesos mediante la operaci´on MPI Allgather, obteniendo en overlaps (array de St MeshRegion con tantos elementos como procesos) el conjunto de intersecciones en todos los 62 CAP´ ITULO 5. DESARROLLOS procesos. El siguiente paso es crear un comunicador MPI para los procesos en los que la intersecci´on calculada no sea nula. Cabe destacar que la comunicaci´on MPI Allgather es necesaria ya que posteriormente se debe llamar aDACreate con los mismos argumentos en todos los procesos y, por tanto, los c´alculos que siguen se deben realizar en todos los procesos implicados. count = 0; for (i = 0; i < mpi_size; i++) if (overlaps[i].totalCells > 0) ranks[count++] = i; ... MPI_Comm_group(PETSC_COMM_WORLD, &worldGroup); MPI_Group_incl(worldGroup, count, ranks, &patchGroup); MPI_Comm_create(PETSC_COMM_WORLD, patchGroup, &patchComm); De esta manera se crea un nuevo grupo conteniendo ´unicamente los procesos implicados (ranks contiene el rango de los procesos con intersecci´on no nula) y un nuevo comunicador para ese grupo, que posteriormente se almacenar´a en la estructura St Patch. Posteriormente, es necesario calcular cu´antas divisiones horizontales y verticales tendr´a el patch y cu´al ser´a el tama˜no de ´estas para que se correspondan con las intersecciones obtenidas: divs = 0; for (i = 0; i < mpi_size; i++) if (overlaps[i].totalCells > 0) { found = PETSC_FALSE; for (j = 0; j < i && !found; j++) if (overlaps[j].totalCells > 0 && overlaps[i].lowerCorner[dir] == overlaps[j]. lowerCorner[dir]) { found = PETSC_TRUE; } if (!found) sizes[divs++] = overlaps[i].nCells[dir]; } siendo dir el eje considerado en cada caso. As´ı, repiti´endose este proceso tres veces, para cada direcci´on se obtiene el n´umero de coordenadas de inicio distintas para todas las intersecciones, junto con el n´umero de celdas a lo largo de dicho eje para cada divisi´on distinta. Por ´ultimo, s´olo queda la creaci´on del objeto DA: if (overlaps[Local.rank].totalCells > 0) { DACreate3d(patchComm, periodicity, DA_STENCIL_BOX, nx, ny, nz, m, n, p, dof, stencil, lx, ly, lz, &da); } else { 5.1. MANTENIMIENTO Y OPTIMIZACI ´ ON DEL C ´ ODIGO 63 da = PETSC_NULL; } donde periodicity,dof ystencil son par´ametros de la creaci´on para la periodicidad, el n´umero de grados de libertad y el n´umero de nodos ghost respectivamente; m,nypson el n´umero de divisiones obtenido anteriormente en divs para cada direcci´on; y lx,ly ylz son los tama˜nos de las divisiones obtenidos anteriormente en sizes para cada direcci´on. De esta manera, se obtiene un DA para el patch en todos aquellos procesos que deben manejar parte del mismo. Adem´as, en todos aquellos procesos donde la intersecci´on del patch con el subdominio local sea nula, el objeto DA ser´a nulo y se podr´a utilizar esta comprobaci´on para, en cualquier otro punto del c´odigo, saber si el patch es local o no. Esta comprobaci´on es de vital importancia ya que las posteriores llamadas a las funciones de tipo DAXXX deber´an realizarse ´unicamente por aquellos procesos cuyo objeto DA no sea nulo. 5.1.4. Refactorizaci´on de las propiedades Como se ha comentado en el apartado anterior, exist´ıa un comportamiento divergente en ciertos casos y que daba resultados notablemente diferentes para casos paralelos con respecto a los resultados secuenciales. Puesto que las propiedades se evaluaban de manera local con nodos vecinos y no exist´ıa comunicaci´on entre ellas, se pens´o que pod´ıa ser una causa de dicho comportamiento. Adem´as, puesto que la evaluaci´on en local sin nodos vecinos pod´ıa solucionar este problema y facilitar otras implementaciones, se opt´o por incluir objetos Vec en la implementaci´on de la estructura St Prop (en lugar de los arrays double* iniciales) y utilizarlos, junto con los DA explicados en el apartado anterior, para realizar la comunicaci´on. As´ı, se incluy´o el atributo values de tipo Vec que se crea de la siguiente manera: if (patch->da != PETSC_NULL) { ... DACreateLocalVector(patch->da, &prop->values); VecSetFromOptions(prop->values); ... } siendo patch (St Patch) la regi´on de la malla donde se define la propiedad y prop (St Prop) la propiedad que se est´a creando. De esta manera, se puede acceder al vector mediante las llamadas de tipo VecXXX o mediante el array obtenido con la llamada a VecGetArray. Cabe destacar que en el caso de acceder a los elementos del vector mediante array, es necesario tener en cuenta que para una propiedad con m´as de un grado de libertad, los elementos en el vector se ordenan de la manera (a0, b0, c0, . . . , a1, b1, c1, . . . , an, bn, cn, . . . ), 64 CAP´ ITULO 5. DESARROLLOS donde la letra indica el grado de libertad y el sub´ındice, el ´ındice de la celda. De esta manera, para acceder al grado de libertad dde la celda ien una propiedad con ngrados de libertad, se debe usar prop->array[i * n + d]. El mismo comportamiento es aplicable a los vectores utilizados para el c´alculo de medias, aunque en principio el acceso a elementos particulares de estos vectores carece de sentido. Por ´ultimo, gracias a la implementaci´on de los patches mediante objetos DA es posible realizar la evaluaci´on de propiedades ´unicamente en las celdas locales de cada proceso (sin nodos ghost) y realizar posteriormente la comunicaci´on de la siguiente manera: DALocalToLocalBegin(prop->Patch->da, prop->values, INSERT_VALUES, prop->values); DALocalToLocalEnd(prop->Patch->da, prop->values, INSERT_VALUES, prop->values); que se encarga de actualizar los nodos ghost del vector values en todos los procesos con los valores correctos, descartando los valores almacenados anteriormente. 5.2. Nuevas funcionalidades Al margen de desarrollos dedicados a la depuraci´on, optimizaci´on y refactorizaci´on del c´odigo existente, tambi´en se han producido desarrollos que dotan a MICSc de nueva funcionalidad. En los siguientes apartados se detallan los m´as importantes. 5.2.1. Introducci´on de la viscosidad Al comienzo del proyecto, la ´unica posibilidad de contemplar la viscosidad propia de fluidos turbulentos (implementada mediante el modelo de Smagorinsky 2.3.2) era introducirla como parte del m´odulo de usuario, que defin´ıa la funci´on de c´alculo de la viscosidad turbulenta a partir de la densidad, la viscosidad propia del modelo LES y la viscosidad especificada como par´ametro de entrada. Tras la refactorizaci´on de MICSc a la estructura de librer´ıa con diferentes casos, la viscosidad se introdujo como una modificaci´on de Gamma mediante la llamada MICScVarSetGamma donde se le especificaba como par´ametro la funci´on de evaluaci´on de la viscosidad. Sin embargo, teniendo en cuenta que entre los principales objetivos del proyecto se encontraba la simulaci´on de flujos turbulentos mediante LES, y que la funci´on de evaluaci´on de la viscosidad es conocida y no var´ıa de un caso a otro, se decidi´o que la viscosidad deb´ıa formar parte del n´ucleo de la librer´ıa para que el usuario pudiera introducirla en su caso de la manera m´as c´omoda posible. 5.2. NUEVAS FUNCIONALIDADES 65 As´ı, se ha tomado una implementaci´on de la funci´on de la viscosidad total proporcionada por A. Cubero y se ha introducido en la librer´ıa de forma que en el momento de la creaci´on de las velocidades, en funci´on de si el usuario especifica un par´ametro de entrada (-les), es posible determinar si se debe introducir el modelo LES y, en su caso, crear las variables especificando la propiedad Gamma con la funci´on de la viscosidad: if (Input.useLES) { CreateProp(DomainPatch, VARI, Input.mu, ViscTotal, 1, CDS, 1, PETSC_FALSE, &gamma); } ... CreateVar("U", CDS, -10, 10.0, linRel, fdtRel, resref, corref, initValues, gamma, &varU); Al margen de valores de relajaci´on y dem´as par´ametros de estas llamadas, es importante destacar que la llamada a la funci´on CreateProp crea la propiedad (St Prop)gamma, que se eval´ua mediante la funci´on ViscTotal. As´ı, en la creaci´on de la variable (CreateVar) se env´ıa como par´ametro la propiedad gamma, que act´ua como coeficiente de difusi´on para la componente xde la velocidad. La implementaci´on de ViscTotal no se incluye en este documento pero sigue el modelo de Smagorinsky y se puede encontrar en el c´odigo fuente y la documentaci´on. Llegados a este punto es importante destacar que la viscosidad afecta a todas las velocidades pero el c´alculo de ´esta no depende de dichas velocidades. Teniendo en cuenta que en un principio cada una de las propiedades asociadas a las variables eran exclusivas de esa variable, y con el fin de optimizar la aplicaci´on tanto en el uso de memoria como en el c´alculo, se ha modificado ligeramente la creaci´on de variables, tal y como se ha mostrado, para que varias variables puedan compartir una misma propiedad de forma que dentro del sistema existir´a una ´unica propiedad para la evaluaci´on de la viscosidad que se evaluar´a una sola vez en cada iteraci´on exterior y se almacenar´a una sola vez donde sea accesible para todas las variables que la necesiten. 5.2.2. C´alculo de medias Como se explica de manera detallada en [22], existen una gran variedad de elementos que pueden llevar a la obtenci´on de resultados instant´aneos distintos para un mismo problema. Estos factores comprenden desde elementos propios de la computaci´on como el n´umero de procesadores, la precisi´on de la m´aquina o los errores de redondeo, hasta factores f´ısicos como las condiciones iniciales o la propia inestabilidad del fen´omeno turbulento. Es por esto 66 CAP´ ITULO 5. DESARROLLOS que para el an´alisis de los resultados proporcionados por la simulaci´on, es interesante calcular ciertos valores estad´ısticos sobre los valores que toman las distintas variables del problema. Estos valores son, tanto las medias de los valores instant´aneos, como las medias de las covarianzas de una variable con respecto de todas las dem´as (incluyendo la propia variable), que no son m´as que las varianza y covarianzas muestrales de las variables. Cabe destacar que los valores sobre los que se han de calcular las medias son siempre aquellos que se han alcanzado en cada iteraci´on de tiempo, al finalizar la convergencia de las iteraciones exteriores de Picard; es decir, no tiene sentido tener en cuenta los valores intermedios tras resolver los sistemas lineales puesto que estos valores no tienen sentido f´ısico mientras no se alcance la convergencia de la iteraci´on de Picard. As´ı pues, para ilustrar los c´alculos que se deben realizar, utilizaremos un ejemplo en dos dimensiones donde las medias a calcular son las siguientes: u: Media de los valores instant´aneos de la velocidad sobre el eje X. Se calcula como 1 k k X i=1 ui(5.1) donde kes el n´umero de pasos de tiempo realizados y uies el valor instant´aneo de la componente X de la velocidad en el paso de tiempo tk. v: Media de los valores instant´aneos de la velocidad sobre el eje Y. Se calcula de manera an´aloga a la componente X sustituyendo uipor vi (5.1). p: Media de los valores instant´aneos de la presi´on. Se calcula de manera an´aloga a las velocidades sustituyendo uipor pi(5.1). u0u0: Varianza muestral de la componente X de la velocidad. Se calcula como: 1 k k X i=1 (ui−u)2(5.2) v0v0: Varianza muestral de la componente Y de la velocidad. Se calcula de manera an´aloga a la componente X sustituyendo uiyupor viyv respectivamente (5.2). p0p0: Varianza muestral de la presi´on. Se calcula de manera an´aloga a las componentes de la velocidad sustituyendo uiyupor piyprespectivamente (5.2). 5.2. NUEVAS FUNCIONALIDADES 67 u0v0: Covarianza muestral de las dos componentes (X e Y) de la velocidad. Se calcula de manera muy similar a la varianza: 1 k k X i=1 (ui−u)(vi−v) (5.3) u0p0: Covarianza muestral de la componente X de la velocidad con la presi´on. Se calcula sustituyendo viyvpor piypen la ecuaci´on anterior (5.3). v0p0: Covarianza muestral de la componente Y de la velocidad con la presi´on. Se calcula sustituyendo uiyupor piypen la ecuaci´on anterior (5.3). Con el fin de calcular las medias, varianzas y covarianzas de las variables, y sabiendo que en cada paso de tiempo el valor de las variables (tras escribirlo a disco en caso necesario) no se mantiene en memoria, es necesario conocer en cualquier momento los sumatorios de los numeradores. Adem´as, puesto que para el c´alculo de la varianza y la covarianza, es necesario conocer la media en cada instante de tiempo, las medias se calcular´an y se almacenar´an en memoria en cada paso de tiempo. Sin embargo, puesto que no en todas las iteraciones temporales se va a escribir a disco y, por tanto, no va a ser necesario conocer todos los valores estad´ısticos, la divisi´on del sumatorio por el n´umero de pasos de tiempo en las ecuaciones de la varianza y covarianza se optimizar´an de forma que ´unicamente se calcular´an cuando se vayan a escribir a fichero. 5.2.3. Modelo de van Driest El modelo de van Driest modifica la constante de Smagorinsky con el fin de amortiguarla en las celdas m´as cercanas a la pared. Para ello, se mostr´o en 2.3.2 que es necesario conocer la distancia a la superficie m´as cercana, junto con la velocidad tangencial a dicha superficie. De esta manera, puesto que este modelo proporciona una clara mejora frente al modelo de Smagorinsky simple, y no plantea ning´un problema adicional, se implement´o este modelo por defecto siempre que se afronte un caso turbulento con LES (-les). Para esta implementaci´on, en primer lugar es necesario calcular y almacenar la distancia a la pared m´as cercana (Geom.wallDistance), junto con un vector para cada una de las celdas (Geom.normalVec). Este vector tendr´a en cuenta la direcci´on desde la propia celda hasta la celda de pared m´as cercana y calcular´a un vector unitario normal a dicha direcci´on, que maximice la componente X y que minimice la componente Y de dicho vector. En la Figura 5.2 se muestra un ejemplo en dos dimensiones. De esta manera, se asume que el fluido siempre fluye de oeste a este, ya que es imposible obtener un vector normal en 3D en la direcci´on y sentido del fluido si no se 74 CAP´ ITULO 5. DESARROLLOS output = DBOpen(localFile, DB_HDF5, DB_APPEND); ... var = firstVar; ... while (var != PETSC_NULL) { WriteMultivar(output, var->name, 1, 0, &var->Array); var = var->next; } if (instantMeans != PETSC_NULL) { CreateAndEnterDir(output, meansPath); var = firstVar; for (i = 0; i < Input.nvar; i++) { PetscSNPrintf(name, 1024, "Mean_%s_", var->name); VecGetArray(instantMeans[i], &array); WriteMultivar(output, name, 1, 0, &array); var = var->next; } ... } DBClose(output); En este caso se obvian tanto la reserva y liberaci´on de memoria, como el caso de la escritura de las medias de las varianzas, que se implementan de manera pr´acticamente id´entica a las medias de los valores instant´aneos. var es una variable de tipo St Var;meansPath es una cadena con la ruta donde se deben escribir los resultados de las medias dentro del fichero SILO; name es una cadena y array es un vector de valores reales; adem´as, CreateAndEnterDir es una funci´on que equivale a mkdir <dir>; cd <dir> en un sistema UNIX. Se puede observar como tanto la escritura de variables como de medias descansa sobre la funci´on WriteMultivar. Esta funci´on recibe una referencia al fichero SILO donde se debe escribir la variable, el nombre, los grados de libertad, una lista de cadenas para los nombres de las subvariables (si el grado de libertad es mayor que uno) y la lista de los valores de las variables. Es importante notar que todos los procesos llaman a esta funci´on, cada uno con los valores locales de la variable. La implementaci´on de WriteMultivar es la siguiente: WriteMultivar(DBfile *output, const char *varName, PetscInt numSubVars, char **subVarNames, PetscReal **values) { ... if (numSubVars == 1) { 5.2. NUEVAS FUNCIONALIDADES 75 DBPutQuadvar1(output, varName, "/Mesh", values[0], Local.meshRegion.nCells, 3, PETSC_NULL, 0, DB_DOUBLE, DB_NODECENT, PETSC_NULL); } else { DBPutQuadvar(output, varName, "/Mesh", numSubVars, subVarNames, values, Local.meshRegion.nCells, 3, PETSC_NULL, 0, DB_DOUBLE, DB_NODECENT, PETSC_NULL); } if (Local.rank == 0) { ... global = DBOpen("output.silo", DB_HDF5, DB_APPEND); CreateAndEnterDir(global, path); for (i = 0; i < mpiSize; i++) { PetscSNPrintf(varNames[i], 1024, "output_%d.silo" ":%s/%s", i, path, varName); varTypes[i] = DB_QUADVAR; } DBPutMultivar(global, varName, mpiSize, varNames, varTypes, PETSC_NULL); DBClose(global); ... } ... } Como se puede observar, el funcionamiento es completamente an´alogo a la implementaci´on de la escritura de la malla. Las diferencias radican en las funciones para escribir las variables (DBPutQuadVar yDBPutMultivar) y en la consideraci´on de que una variable pueda estar formada por distintas subvariables (numSubVars >1), como es el caso de la velocidad que, por simplicidad, no se ha detallado al completo. Por ´ultimo, destacar que el caso de la escritura de propiedades se basa tambi´en en llamadas a la funci´on WriteMultivar y es pr´acticamente id´entico al caso de la escritura de variables, por lo que no se incluye su implementaci´on en este documento. 5.2.6. Configure /Makefile La posibilidad de escribir los resultados en formato SILO implica que se debe tener instalado tanto SILO con soporte para HDF5 y se debe de poder enlazar MICSc con estas librer´ıas en tiempo de compilaci´on con el fin de dar la funcionalidad requerida. Sin embargo, existen ciertos casos en los que esta 76 CAP´ ITULO 5. DESARROLLOS funcionalidad puede no ser posible o incluso contraproducente. Por ejemplo, es posible que se est´e intentando ejecutar en una m´aquina donde no se tengan permisos de administrador o no se tenga acceso a Internet para obtener SILO y HDF5; o bien simplemente se quiere realizar una ejecuci´on para comprobar la correcci´on del c´odigo o medir tiempos de cara a calcular aceleraciones y eficiencia como los presentados en 6.2.2. En estos casos resulta m´as adecuado poder compilar MICSc sin soporte para SILO. As´ı, con el fin de poder realizar distintas configuraciones de MICSc, y adecu´andose a las convenciones GNU [25], se gener´o un script de configuraci´on configure y un fichero makefile con la regla make para la compilaci´on de MICSc. El script configure se encuentra escrito en Python y para su implementaci´on, al igual que para la implementaci´on de ciertas secciones del fichero makefile, se utilizaron como documentaci´on y ejemplo las implementaciones de SLEPc [26], un proyecto llevado a cabo por compa˜neros del mismo grupo de investigaci´on (GRyCAP [27]). As´ı pues, ´unicamente comentaremos las opciones que proporciona este script: --with-silo:Configura la instalaci´on buscando las librer´ıas y cabeceras de SILO en los directorios por defecto (i.e. /usr/local,/opt, ...). --with-silo-dir=<silo dir>:Configura la instalaci´on de MICSc buscando las librer´ıas de SILO en <silo dir>/lib y las cabeceras en <silo dir>/include. --with-silo-flags=<silo flags>:Configura la instalaci´on especificando exactamente los par´ametros que se deben utilizar en el comando de compilaci´on (por ejemplo, -lsiloh5 -lhdf5 -L<dir> -I<dir> ...). Cabe destacar que en el caso de especificar el directorio de SILO, las librer´ıas de HDF5 deber´an estar en los directorios por defecto. En caso de no estar en dichos directorio, la configuraci´on fallar´a. Esto tiene como soluci´on incluir los argumentos de la compilaci´on mediante --with-silo-flags. Sin embargo, aunque efectiva, esta soluci´on es algo inc´omoda y admite ciertas mejoras menores que se contemplan en la hoja de ruta del proyecto pero que no merecen una menci´on especial en la secci´on de trabajo futuro (7.2). Una vez se ha realizado la configuraci´on, se debe definir la variable de entorno MICSC DIR (si no se ha definido antes) que contenga la ruta ra´ız de MICSc. Adem´as, se har´a uso de la variable de entorno PETSC ARCH para crear un nuevo directorio dentro de $MICSC DIR que contendr´a las librer´ıas, cabeceras y ficheros de configuraci´on espec´ıficos. Posteriormente, mediante la orden make se compilar´a la librer´ıa. Llegados a este punto, es posible compilar todos los casos proporcionados con MICSc mediante la regla make examples. 5.2. NUEVAS FUNCIONALIDADES 77 Por ´ultimo, el fichero makefile tambi´en proporciona diversas reglas independientes que nos permiten, entre otras cosas, crear un conversor de SILO a Tecplot, de binario a Tecplot o una aplicaci´on basada en PETSc para la comparaci´on de normas de vectores en distintos ficheros binarios. La implementaci´on de esto se basa principalmente en las pr´acticas que utiliza PETSc para su configuraci´on y compilaci´on y se puede encontrar en la documentaci´on y el c´odigo fuente. De esta manera, y teniendo en cuenta que la ´unica herramienta de compilaci´on al principio del proyecto era un fichero makefile para la compilaci´on de todo el proyecto (con los problemas planteados en 5.1.1) con una sola regla main, estas implementaciones tambi´en constituyen uno de los avances en usabilidad m´as importantes de este proyecto. 78 CAP´ ITULO 5. DESARROLLOS Cap´ıtulo 6 Experimentos y resultados En este cap´ıtulo se detallan los experimentos realizados con el c´odigo para su presentaci´on en el congreso ECT2010 celebrado en Valencia. En ellos se analizan tanto los resultados de cara a validarlos frente a otros resultados conocidos, como las medidas temporales con el fin de obtener una estimaci´o´n del comportamiento y eficiencia paralelos. Cabe destacar que existen otros casos de validaci´on realizados anteriormente para el c´odigo en el contexto de la tesis doctoral de A. Cubero que pueden encontrarse en [14]. 6.1. Validaci´on LES: caso del canal peri´odico Para la validaci´on del comportamiento del c´odigo en la resoluci´on de flujos turbulentos mediante LES se utiliz´o el caso del canal peri´odico. Este caso consiste en un flujo turbulento a trav´es de un canal caracterizado por tener dos condiciones de contorno de tipo pared en las regiones norte y sur; es decir, a lo largo del plano XZ para los valores m´ınimo y m´aximo de Y en el dominio. Adem´as, se caracteriza por la periodicidad tanto en el eje X como en el eje Y. En nuestro caso concreto, la malla utilizada consta de 121 ×121 ×81 nodos a lo largo de las direcciones x, y, z respectivamente. Por otro lado, cabe destacar que se ha utilizado una malla con diferente tama˜no en los vol´umenes finitos a lo largo de la direcci´on normal al flujo (eje Y), siguiendo una determinada funci´on hiperb´olica. El resultado es una malla representando el dominio [0, 0, 0] - [6.45, 2, 3.24] como se muestra en la Figura 6.1. Para la validaci´on, en primer lugar se realiz´o la ejecuci´on del caso con la implementaci´on simple del modelo de Smagorinsky para la turbulencia, como se detalla en 2.3.2. Adem´as, con el fin de contrastar y validar los resultados obtenidos, estos se compararon con dos simulaciones distintas: una realizada por Moser [28] siguiendo la aproximaci´on DNS y otra realizada por Piomelli mediante LES. Por ´ultimo, con el fin de obtener un flujo con el n´umero de Reynolds adecuado y poder comparar el flujo frente a estas otras simulaciones, se fij´o un diferencial de presi´on dp dx cons79 80 CAP´ ITULO 6. EXPERIMENTOS Y RESULTADOS Figura 6.1: Malla 121×121×81 para el caso del canal peri´odico. tante e igual a 1. Los resultados obtenidos se muestran en la Figura 6.2. En estos resultados se puede observar como los perfiles de variables promediadas salen cualitativamente bien, en el sentido de que las curvas son suaves, sim´etricas y con los picos localizados en las regiones donde se espera. Sin embargo, las comparaciones con los resultados de Moser y Piomelli no obtienen el resultado deseado, debido a las grandes diferencias en magnitud entre simulaciones. En este momento se pens´o que las diferencias eran debidas al modelo LES utilizado, y que la implementaci´on del modelo de van Driest 2.3.2 que amortigua el efecto de la turbulencia en las celdas cercanas a las paredes podr´ıa dar mayor precisi´on a los resultados. As´ı pues, tras realizar la implementaci´on del modelo de van Driest, y siendo ya presentado en el congreso ECT2010 celebrado en Valencia [17], se obtuvieron los resultados que se muestran en la Figura 6.3. En estos resultados se puede observar como la tendencia de los resultados, comparados esta vez ´unicamente con el DNS de Moser, sigue siendo satisfactoria, pero adem´as la magnitud de los resultados se acerca hasta alcanzar unos valores aceptables. Sin embargo, los resultados de la simulaci´on muestran como el flujo se desarrolla hasta alcanzar el estado estad´ısticamente estacionario con un n´umero de Reynolds de 224, mientras que los resultados de Moser se aplican a un flujo con un n´umero de Reynolds de 180. En la pr´actica, esto resulta en que se est´a validando una simulaci´on frente a otra muy similar pero no id´entica. De esta manera, las diferencias menores en los resultados pueden corregirse (o por el contrario, acrecentarse) simulando con el mismo n´umero 6.2. RENDIMIENTO DEL C ´ ODIGO PARALELO 81 0 5 10 15 20 25 30 1 10 100 u+ y+ - 1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 uv+ y LES Piomelli MICSc DNS Moser LES Piomelli MICSc DNS Moser Figura 6.2: Resultados para la simulaci´on del caso del canal peri´odico con el modelo de Smagorinsky. 0 5 10 15 20 25 1 10 100 u+ y+ DNS Moser MICSc -1 -0.5 0 0.5 1 0 0.5 1 1.5 2 uv+ y DNS Moser MICSc Figura 6.3: Resultados para la simulaci´on del caso del canal peri´odico con el modelo de van Driest. de Reynolds. As´ı, para obtener el n´umero de Reynolds deseado, el diferencial de presi´on fijado deber´ıa de ajustarse en cada paso de tiempo para que de esta manera no tuviera que variar la velocidad media del fluido, y as´ı mantener el n´umero de Reynolds constante e igual a 180. 6.2. Rendimiento del c´odigo paralelo Para cuantificar el rendimiento del c´odigo en el caso del canal peri´odico, explicado en la secci´on anterior, en esta secci´on se presentan ciertas medidas temporales, principalmente para la discusi´on de dicho rendimiento en t´erminos de eficiencia paralela. En cuanto a la configuraci´on utilizada para las simulaciones, se realizaron 30 pasos temporales con un total de 360 resoluciones de sistemas lineales en el total de cada una de las pruebas. Las pruebas se realizaron en 82 CAP´ ITULO 6. EXPERIMENTOS Y RESULTADOS KSP PC Setup Solve Total Its gmres jacobi 0.019 0.692 0.711 12.84 gmres bjacobi 0.129 0.322 0.451 3.83 bcgs jacobi 0.019 0.570 0.589 7.25 bcgs bjacobi 0.127 0.332 0.459 2 bcgsl jacobi 0.020 0.812 0.832 9.33 bcgsl bjacobi 0.129 0.556 0.685 3.67 tfqmr jacobi 0.019 1.052 1.071 11.86 tfqmr bjacobi 0.130 0.353 0.483 2.03 Tabla 6.1: Comparaci´on de solvers lineales (KSP) y precondicionadores (PC). CaesarAugusta, un supercomputador con 256 nodos JS20, cada uno de ellos con dos procesadores PowerPC 970FX de 64 bits funcionando a 2.2 GHz, interconectados con una red Myrinet de baja latencia. 6.2.1. Comparaci´on de solvers lineales y precondicionadores La primera prueba se centra en la comparaci´on en cuanto a tiempo y n´umero de iteraciones en la resoluci´on de los sistemas lineales que componen el problema. Concretamente, se considera el tiempo medio de setup, de resoluci´on y el total (suma de ambos), todo ello en valores medios por resoluci´on de sistema lineal. Por otro lado, el n´umero de iteraciones se muestra de la misma manera, promediado por cada resoluci´on de sistema lineal, que como se ha comentado anteriormente suman un total de 360 resoluciones. En cuanto a los solvers lineales utilizados para la comparaci´on, estos son GMRES (gmres), BiCGStab (bcgs) y BiCGStab(`) con `= 2 (bcgsl), mientras que los precondicionadores utilizados son Jacobi (jacobi) y Jacobi por bloques o Block Jacobi (bjacobi). En el caso de Jacobi por bloques se utilizan pbloques, siendo pel n´umero de procesadores utilizados en la simulaci´on y utilizando una factorizaci´on incompleta LU o ILU por cada uno de los bloques. En el caso concreto de esta prueba, el n´umero de procesadores utilizados para todas las ejecuciones fue de 32. Se pueden encontrar m´as detalles sobre estos m´etodos y precondicionadores en [29]. Los resultados obtenidos se muestran en la Tabla 6.1. En esta tabla se puede observar en primer lugar que el uso de Block Jacobi reduce de manera notable el n´umero de iteraciones necesarias para alcanzar la convergencia, mientras que por otro lado aumenta sensiblemente el tiempo de setup. Esto es debido a que, por un lado el precondicionador de Jacobi se obtiene de forma inmediata a partir de la diagonal de la matriz Adel sistema Ax =b, mientras que el precondicionador por bloques requiere de la computaci´on de la factorizaci´on LU para cada uno de los bloques diagonales que, aunque incompleta, obviamente requiere un mayor tiempo de computaci´on. 6.2. RENDIMIENTO DEL C ´ ODIGO PARALELO 83 Por otro lado, el tiempo por iteraci´on con el uso de Jacobi por bloques es mayor que para Jacobi simple, pero siendo el n´umero de iteraciones menor, en la medici´on del tiempo global, el resultado es siempre menor con el uso de Jacobi por bloques. Por ´ultimo, en cuanto a la comparaci´on de solvers lineales, los mejores resultados los presentan tanto GMRES como BiCGStab, con pr´acticamente los mismos resultados. Es por esto que para la evaluaci´on de aceleraci´on y eficiencia del siguiente apartado se consideran las combinaciones GMRES + Block Jacobi y BiCGStab + Block Jacobi, junto con la combinaci´on que exist´ıa por defecto cuando se adopt´o el c´odigo, GMRES + Jacobi simple. 6.2.2. Evaluaci´on de aceleraci´on y eficiencia paralela Como se comenta en el apartado anterior, en esta prueba se realiz´o la ejecuci´on del caso del canal peri´odico con distintas configuraciones de solvers lineales y precondicionadores. Concretamente GMRES + Block Jacobi, BiCGStab + Block Jacobi y GMRES + Jacobi simple. Para cada una de estas pruebas se realizaron varias pruebas con distinto n´umero de procesadores p, en concreto p={1,2,4,8,16,32,64}. Adem´as, para cada ejecuci´on se muestran dos secciones: en la primera de ellas se muestra el tiempo total de la ejecuci´on junto con la aceleraci´on (Sp) y eficiencia (Ep), mientras que en la segunda se considera ´unicamente el tiempo medio de resoluci´on de sistemas lineales (dividido en setup y resoluci´on), tambi´en acompa˜nado de las correspondientes aceleraci´on (Sp) y eficiencia (Ep). Los resultados se muestran en la Tabla 6.2. De esta tabla se puede extraer una conclusi´on global a partir de los resultados de eficiencia para los tiempos globales, que incluyen tanto la resoluci´on de sistemas como la evaluaci´on previa de propiedades y el ensamblado del sistema. Estos resultados, muestran una tendencia a descender conforme se aumenta el n´umero de procesadores, tal y como se esperaba; pero por otro lado, el valor de la eficiencia se muestra por encima del 60 % para todos los casos, incluso para el caso de 64 procesadores. Es por ello que los valores de aceleraci´on y eficiencia para el tiempo total se pueden considerar razonablemente satisfactorios. Sin embargo, para el caso de los solvers lineales, los resultados son ligeramente peores. En concreto, BiCGStab muestra una eficiencia muy a la par con la eficiencia global de la ejecuci´on, mientras que el peor resultado se muestra para la combinaci´on de GMRES + Jacobi, estando siempre la aceleraci´on y eficiencia siempre por debajo de los resultados del tiempo total, muy posiblemente debido al mayor n´umero de iteraciones efectuadas. Por ´ultimo, se puede observar que el comportamiento del precondicionador de Block Jacobi es satisfactorio, a pesar de que la efectividad disminuye conforme aumenta el n´umero de procesadores, ya que el tama˜no del bloque es cada vez menor. 90 AP´ ENDICE A. MANUAL DE USUARIO DE MICSC --with-silo-flags=<silo flags>:Configura la instalaci´on especificando exactamente los par´ametros que se deben utilizar en el comando de compilaci´on (por ejemplo, -lsiloh5 -lhdf5 -L<dir> -I<dir> ...). De esta manera, es posible realizar la configuraci´on con soporte para SILO: $ ./configure --with-silo-dir=<silo_dir> o bien realizar una instalaci´on sin soporte para SILO, con PETSc como ´unica dependencia. $ ./configure Llegados a este punto, es imprescindible definir la variable de entorno MICSC DIR (en caso de no haberse definido ya anteriormente) para poder proceder con la compilaci´on. $ MICSC_DIR=$PWD; export MICSC_DIR $ make Esta compilaci´on crear´a un directorio $PETSC ARCH dentro del directorio de MICSc ($MICSC DIR) que contendr´a la librer´ıa compilada as´ı como cualquier otro fichero espec´ıfico de la configuraci´on realizada. En este punto se puede comprobar la correcci´on de la instalaci´on mediante $ make test Es importante destacar que estos tests hacen uso MPI, por lo que el testeo en m´aquinas con sistemas de colas para este tipo de ejecuciones puede ser problem´atico. Adem´as, cabe destacar que es posible compilar los casos b´asicos proporcionados con la distribuci´on de MICSc con $ make examples Los binarios generados, junto con el resto de ficheros necesarios para la ejecuci´on de estos casos se pueden encontrar en $MICSC DIR/cases, en un subdirectorio distinto para cada uno de los casos. Ejecuci´on MICSc proporciona diferentes casos de uso sobre los que se pueden realizar diferentes pruebas. Para ello, accedemos al directorio de un caso como, por ejemplo, el del canal peri´odico, lo compilamos (si no lo hemos compilado ya antes con make examples) y lo ejecutamos. 91 $ cd $MICSC_DIR/cases/periodchannel $ make $ mpirun -np <p> ./main <args> Las opciones que admiten los casos generados con MICSc se pueden obtener con el argumento -help. Al margen de las opciones proporcionadas por PETSc para vectores, matrices, solvers, precondicionadores y dem´as, las opciones de MICSc son las siguientes (<n> determina un entero, run real, eun enumerado y fun fichero): Opciones de inicializaci´on -grid file <f> : Fichero de malla. -rod file <f> : Fichero con los par´ametros por defecto del caso. Opciones de monitorizaci´on -log level <e> : Nivel de monitorizaci´on {0, 1, 2} -out ievery <n> : N´umero de iteraciones exteriores a realizar para escribir resultados. -1 para no escribir resultados intermedios. -out idtevery <n> : N´umero de pasos temporales a realizar para escribir resultados. -1 para no escribir resultados intermedios. -start instant means <n> : Paso de tiempo en el que se debe empezar el c´alculo de medias de valores instant´aneos. -1 para no realizar dicho c´alculo. -start variance means <n> : Paso de tiempo en el que se debe empezar el c´alculo de medias de las varianzas. -1 para no realizar dicho c´alculo. -inifl : Determina si se debe escribir el resultado inicial. -finfl : Determina si se debe escribir el resultado final. -print props : Determina si se deben escribir resultados de propiedades. Actualmente la ´unica propiedad que admite la escritura de resultados es la viscosidad LES. Opciones de resoluci´on del sistema de ecuaciones -max outer its <n> : N´umero m´aximo de iteraciones exteriores. -min outer its <n> : N´umero m´ınimo de iteraciones exteriores. -res max <r> : Norma m´axima del residuo normalizado. -res max <r> : Norma m´axima del vector correcci´on normalizado. -k restart : Determina si se debe reiniciar una ejecuci´on anterior a partir de un fichero de resultados binario restart.bin. 92 AP´ ENDICE A. MANUAL DE USUARIO DE MICSC -init it <n> : N´umero de la primera iteraci´on. -k wrap <e> : Periodicidad del dominio del problema {NONPERIODIC, XPERIODIC,YPERIODIC,XYPERIODIC,XYZPERIODIC,XZPERIODIC, YZPERIODIC,ZPERIODIC,XYZGHOSTED} -n ghost <n> : N´umero de nodos vecinos. -gravity <r> : Fuerza de la gravedad. 0 para ignorarla. -n order trans <e> : Esquema temporal {steady,eul1,eul2,adm2, adm3 } -deltat <r> : Paso de tiempo. -tinit <r> : Instante de tiempo inicial. -tfinal <r> : Instante de tiempo final. Opciones de variables del problema -uref <r> : Valor constante para la entrada de flujo en la direcci´on X. -mu <r> : Valor de µ(viscosidad din´amica). -cs0 <r> : Constante de Smagorinsky para el modelo LES. -dpdx <r> : Valor del flujo de masa ∂p/∂x -tref <r> : Valor de referencia para la temperatura. -rho <r> : Valor de ρ(densidad). -ndim <n> : N´umero de dimensiones. -ncoupvar <n> : N´umero de variables acopladas. -cont eq <e> : Ecuaci´on de continuidad {poisson,simple,simplec, simpler } -inlet <e> : Funci´on para el flujo de entrada {u,f,rho } -les : Determina si se debe utilizar el modelo LES. -kemod : Determina si se debe utilizar el modelo K. -var t : Determina si se debe considerar la temperatura como una variable a resolver. Implementaci´on de nuevos casos Al margen de los casos base proporcionados, MICSc permite la creaci´on de nuevos casos. Para ello, en primer lugar se debe crear un nuevo directorio para el caso $ cd $MICSC_DIR/cases $ mkdir <caso> 93 Este directorio generalmente suele contener (al margen de otros ficheros que puedan ser necesarios para el caso particular) cuatro ficheros principales: grd.inp: Contiene la malla y las condiciones de contorno tal y como se especifica en el siguiente apartado. rod.inp: Contiene las opciones por defecto del caso. Se proporciona por comodidad para que no sea necesario especificarlas siempre por l´ınea de comandos. Cabe destacar que las opciones que luego se puedan especificar en el momento de la ejecuci´on sobreescriben las de este fichero. makefile: Fichero para la compilaci´on del caso. Se puede copiar de cualquier otro caso. main.c: C´odigo fuente para la ejecuci´on del caso. Aparte de las funciones propias del caso, consta de una funci´on main en la que se puede hacer uso de las siguientes llamadas a la librer´ıa de MICSc: MICScFinalize(): Libera toda la memoria utilizada por la aplicaci´on y llama a la rutina de finalizaci´on de los objetos de PETSc. En el c´odigo, debe ser la ´ultima llamada. MICScInitialize(PetscInt argc, char** argv[]): Se encarga de realizar la inicializaci´on de PETSc, as´ı como de leer los ficheros de entrada y reservar e inicializar todos los valores y estructuras de datos comunes a todos los casos. Debe ser la primera llamada del programa y los par´ametros de entrada deben de ser punteros a los par´ametros de entrada de la funci´on main. MICScSetProp(St_Prop *prop, St_Patch *patch, ktypeProp propType, PetscReal cons, PetscErrorCode(funct*)(St_Prop *prop), ktypeIp interpType, PetscReal linRelax, PetscTruth calcMeans): Modifica los valores por defecto de una propiedad. El primer par´ametro es el puntero a la propiedad a modificar y el resto son la regi´on de la malla sobre la que se aplica (patch), el tipo de propiedad (propType), el valor de la propiedad si es constante (cons), la funci´on de evaluaci´on si es variable (funct), el tipo de interpolaci´on (interpType) y la relajaci´on lineal (linRelax). Adem´as, es posible especificar si se desean calcular medias para esta propiedad mediante calcMeans. MICScSolve(): 94 AP´ ENDICE A. MANUAL DE USUARIO DE MICSC Resuelve el sistema de ecuaciones. Es el n´ucleo principal de MICSc y la llamada que consume pr´acticamente todo el tiempo de ejecuci´on de la simulaci´on. MICScVarSetGamma(PetscInt index, ktypeProp gammaType, ktypeIp gammaInterpType, PetscReal gammaConst, PetscErrorCode (gammaFunct*) (St_Prop *prop)): Sustituye la propiedad del coeficiente de difusi´on por una nueva con los valores especificados como par´ametro. El primer par´ametro es el ´ındice de la variable; es posible utilizar las macros U,V,W,P, etc para este par´ametro. Esto es aplicable para todas las funciones de nombre MICScVarSet*. MICScVarSetInitValues(PetscInt index, ktypeProp initValuesType, ktypeIp initValuesInterpType, PetscReal initValuesConst, PetscErrorCode (initValuesFunct*) (St_Prop *prop)): Se define de manera an´aloga a MICScVarSetGamma para los valores iniciales de la variable en cada celda. MICScVarSetInterpolation(PetscInt index, ktypeIp interpType): Modifica el tipo de interpolaci´on de la variable. MICScVarSetLimits(PetscInt index, PetscReal minValue, PetscReal maxValue): Modifica los valores m´ınimo y m´aximo de la variable. MICScVarSetReference(PetscInt index, PetscReal resRef, PetscReal corrRef): Modifica los valores de referencia de la variable. MICScVarSetRelaxation(PetscInt index, PetscReal linRelax, PetscReal fdtRelax): Modifica los valores de relajaci´on de la variable. MICScSetInletFunction(PetscErrorCode (*function)(PetscInt*, PetscReal*)): Establece la funci´on que se debe utilizar para calcular el flujo de entrada, en caso de que se defina tal condici´on de contorno. Formato de entrada de la malla Cada caso debe disponer de un fichero que defina tanto la malla como las condiciones de contorno. Este fichero se estructura de la siguiente manera: 95 <{0,1}> <num_divisiones_eje_X> <num_divisiones_eje_Y> <num_divisiones_eje_Z> <lista_coords_eje_X> <lista_coords_eje_Y> <lista_coords_eje_Z> <lista_condiciones_de_contorno> El primer valor representa si la malla es rectangular (0) o cil´ındrica (1). Aunque actualmente s´olo se soportan mallas rectangulares, tambi´en es posible que en un futuro se implementen mallas cil´ındricas. Las siguientes tres l´ıneas contienen un ´unico n´umero entero por l´ınea, que define el n´umero de divisiones para cada uno de los ejes. Posteriormente deben aparecer tres l´ıneas con una lista de coordenadas para cada uno de los ejes. El tama˜no de la lista se deber´a corresponder con el n´umero definido en las l´ıneas anteriores y deber´an estar formados por n´umeros reales separados por espacios. En este punto es importante destacar que se asume la coordenada [0,0,0] como origen de la malla y no es necesario incluirla en este fichero. Por ´ultimo, aparece una lista de condiciones de contorno donde cada una de ´estas ocupa una l´ınea distinta con el siguiente formato: <x> <nx> <y> <ny> <z> <nz> <vecino> <tipo> <porosidad> Los primeros 6 par´ametros son de tipo entero y especifican la regi´on de la malla sobre la que se define la condici´on de contorno ([x, y, z]−[x+ nx, y +ny, z +nz]); el siguiente par´ametro especifica sobre qu´e cara de las celdas se define la condici´on ({north, south, east, west, high, low, none}); el tipo define la condici´on de contorno ({wall, moving wall, inlet, outflow, outpress newman, outpress extrap, symmetry, fix}); y la porosidad es un n´umero real en el rango [-1, 1]. Formato de salida SILO Por ´ultimo, en caso de haber configurado la instalaci´on de MICSc con SILO, los ficheros producidos podr´an abrirse en diversos visores, tales como VisIt. Sin embargo, es necesario conocer la estructura jer´arquica que presentan los ficheros SILO generados por MICSc. Dicha estructura se puede encontrar en la Figura A.1. En esta figura, los elementos marcados en negrita son hojas (variables y malla), mientras que los elementos en cursiva son directorios. As´ı, los resultados se dividen en iniciales (InitialResults), intermedios (IntermediateResults), finales (FinalResults) y de paso temporal (TimeStepResults). Cabe destacar que los resultados de paso temporal 96 AP´ ENDICE A. MANUAL DE USUARIO DE MICSC / Mesh InitialResults InstantValues Velocity P ... FinalResults InstantValues Velocity P ... IntermediateResults field 10 InstantValues Velocity P ... field 20 InstantValues Velocity P ... ... TimeStepResults time 0 01 InstantValues Velocity P Means Mean U Mean V Mean W Mean P ... VarianceMean Velocity VarianceMean P VarianceMean U U VarianceMean U V ... time 0 02 ... ... Figura A.1: Estructura jer´arquica de un fichero SILO generado por MICSc. 97 se escriben al obtener la convergencia (o alcanzar el n´umero m´aximo de iteraciones), mientras que los resultados intermedios se refieren a aquellos resultados producidos en una determinada iteraci´on, a mitad de un paso temporal. Dentro de cada uno de estos directorios habr´a subdirectorios para el paso de tiempo (time <tiempo>) o la iteraci´on (field <it>) cuando corresponda. Por ´ultimo, los directorios m´as alejados de la ra´ız contendr´an valores instant´aneos para las variables (InstantValues) o las medias cuando corresponda (Means). Es importante notar que para las medias Mean <var> indica la media de los valores instant´aneos de la variable var, VarianceMean <var> indica la media de la varianza de la variable var y VarianceMean <var1> <var2> indica la media de la covarianza entre las variables var1 yvar2. 98 AP´ ENDICE A. MANUAL DE USUARIO DE MICSC Bibliograf´ıa [1] S.V. Patankar and D.B. Spalding. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. International Journal of Heat and Mass Transfer, 15:1787–1806, 1972. [2] M. Snir, S. Otto, S. Huss-Lederman, D. Walker, and J. Dongarra. MPI: The Complete Reference. MIT Press, Cambridge, MA, USA, 1995. [3] BLAS - Basic Linear Algebra Subprograms, 2010. http://www. netlib.org/blas. [4] LAPACK - Linear Algebra PACKage, 2010. http://www.netlib.org/ lapack. [5] Satish Balay, Kris Buschelman, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Barry F. Smith, and Hong Zhang. PETSc Web page, 2010. http://www.mcs.anl.gov/ petsc. [6] S. Pope. Turbulent Flows. Cambridge University Press, Cambridge, UK, 2000. [7] J.H Ferziger and M. Peri´c. Computational Methods for Fluid Dynamics. Springer, 3rd edition, 2001. [8] J. Bardina, J.H. Ferziger, and W.C. Reynolds. Improved subgrid scale models for large eddy simulation. AIAA paper, 80-1357, 1980. [9] E.R. van Driest. On turbulent flow near a wall. AIAA journal, 23(11):1007–1010, 1036, 1956. [10] M.T. Heath. Scientific Computing. An Introductory Survey. Mc Graw Hill, 2nd edition, 2000. [11] V. Faber and T. Manteuffel. Necessary and sufficient conditions for the existence of a conjugate gradient method. SIAM Journal on Numerical Analysis, 21(2):352–362, 1984. 99