PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

73

Tercera Parte: CÓDIGO DEL PROYECTO

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

74

I. PROBLEMA FÍSICO: ECUACIÓN DEL CALOR

1. DESCRIPCIÓN DEL PROBLEMA FÍSICO:

La ecuación del calor es una importante ecuación diferencial en

derivadas parciales que describe la distribución del calor (o

variaciones de la temperatura) en una región a lo largo del

transcurso del tiempo. Para el caso de una función de tres variables

en el espacio (x,y,z) y la variable temporal t, la ecuación del calor

es:

donde k es una constante.

La ecuación del calor es de una importancia fundamental en numerosos

y diversos campos de la ciencia.

La resolución de la ecuación del calor mediante métodos numéricos es

un buen ejemplo de la aplicación de esta disciplina matemática a los

problemas reales de la Física, que sin tener una solución analítica,

o teniéndola, pero siendo muy complicada, hacen necesario que se

recurra a una resolución numérica.

2. ANÁLISIS TEÓRICO DEL PROBLEMA:

El problema queda descrito mediante el presente enunciado:

, en Ω

,

.

Con los siguientes datos concretos:

- f=0.

- Ω es el cuadrado [0,1]x[0,1]

- Se identifican los subdominios frontera como sigue:

-I es el lado superior.

-II es el lado derecho.

-III es el lado inferior.

-IV es el lado izquierdo.

-La solución exacta es:

,

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

75

- Se definen como lados Dirichlet el III y el IV y las

condiciones de contorno en esos lados cumplen:

- ,

-

- Se definen como lados Neumann los lados el I y el II y las

condiciones de contorno en esos lados cumplen:

-

-

En la Figura 1, se representa el sistema:

Fig. 1: Representación del sistema objeto de resolución en el dominio

elegido y las condiciones de contorno.

Sol uci ón exact a:

0 1

Lado I I I condi ci ones de

cont or no Di r i chl et

Lado I I condi ci ones de

cont or no Neumann

Lado I V condi ci ones de

cont or no Di r i chl et

Lado I condi ci ones de cont or no

Neumann

1

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

76

En la figura 2 se representa la solución exacta del sistema.

Fig. 2: Representación de la solución exacta.

Se observa que la función solución varía entre -1 y 1. Por la

definición y condiciones de contorno, obtenemos que los lados III y

IV son constantes, mientras que el lado I es responsable del

“enfriamiento” debido al gradiente negativo, y el II es responsable

del “calentamiento“ debido a su gradiente positivo.

Mediante un estudio numérico inicial, detallado en las tablas 1 y 2,

se intenta inferir algunas características, especialmente relevantes

para la sensibilidad de la función que puedan ser útiles en las fases

siguientes del proyecto.

En la tabla 1 se han elegido valores dispares para determinar la

magnitud de la solución con respecto a las coordenadas y por tanto

acotar el error absoluto tolerable, incluyendo valores cercanos al 0.

Como se observa, la mayor sensibilidad se detecta en dicha región y

en la región de la recta diagonal.

En la tabla 2 se han elegido valores concentrados alrededor de los

valores centrales de coordenadas para observar las pautas de cambio

en la salida cuando los nodos están muy cercanos.

Se observa que la solución es aproximadamente del mismo orden de

magnitud que la diferencia entre las coordenadas. Las diferencias en

coordenadas entre dos puntos contiguos determinarán el mínimo y

vendrán fijadas por el mallado y por tanto condicionarán la cota de

error exigible al Solver.

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

77

Como conclusión interesante es que la solución es más sensible al

error en las cercanías de la diagonal x=y y en las cercanías del 0.

Los valores de las coordenadas en la primera región son muy

parecidos y la solución por tanto es pequeña en magnitud; anulándose

en la propia diagonal En la segunda región, las coordenadas son muy

pequeñas y por tanto lo es también la solución. Por ello, una

representación satisfactoria deberá tener suficiente precisión en

ambas regiones.

En las demás regiones se puede tolerar un error superior ya que la

diferencia entre los valores de las coordenadas implicadas es

suficientemente grande para compensarlo.

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

78

Tabla 1: Valores numéricos de la solución en el dominio [0;1]x[0;1].

X/Y 0,0001 0,001 0,01 0,1 0,2 0,3 0,4 0,5 1 0 0,000000 0,000001 0,000100 0,010000 0,040000 0,090000 0,160000 0,250000 1,000000 0,001 -0,000001 0,000000 0,000099 0,009999 0,039999 0,089999 0,159999 0,249999 0,999999 0,01 -0,000100 -0,000099 0,000000 0,009900 0,039900 0,089900 0,159900 0,249900 0,999900 0,1 -0,010000 -0,009999 -0,009900 0,000000 0,030000 0,080000 0,150000 0,240000 0,990000 0,4155 -0,172640 -0,172639 -0,172540 -0,162640 -0,132640 -0,082640 -0,012640 0,077360 0,827360 0,416 -0,173056 -0,173055 -0,172956 -0,163056 -0,133056 -0,083056 -0,013056 0,076944 0,826944 0,417 -0,173889 -0,173888 -0,173789 -0,163889 -0,133889 -0,083889 -0,013889 0,076111 0,826111 0,418 -0,174724 -0,174723 -0,174624 -0,164724 -0,134724 -0,084724 -0,014724 0,075276 0,825276 0,419 -0,175561 -0,175560 -0,175461 -0,165561 -0,135561 -0,085561 -0,015561 0,074439 0,824439 0,4199 -0,176316 -0,176315 -0,176216 -0,166316 -0,136316 -0,086316 -0,016316 0,073684 0,823684 0,42 -0,176400 -0,176399 -0,176300 -0,166400 -0,136400 -0,086400 -0,016400 0,073600 0,823600 0,44 -0,193600 -0,193599 -0,193500 -0,183600 -0,153600 -0,103600 -0,033600 0,056400 0,806400 0,46 -0,211600 -0,211599 -0,211500 -0,201600 -0,171600 -0,121600 -0,051600 0,038400 0,788400 0,48 -0,230400 -0,230399 -0,230300 -0,220400 -0,190400 -0,140400 -0,070400 0,019600 0,769600 0,49 -0,240100 -0,240099 -0,240000 -0,230100 -0,200100 -0,150100 -0,080100 0,009900 0,759900 0,499 -0,249001 -0,249000 -0,248901 -0,239001 -0,209001 -0,159001 -0,089001 0,000999 0,750999 0,4999 -0,249900 -0,249899 -0,249800 -0,239900 -0,209900 -0,159900 -0,089900 0,000100 0,750100 0,5001 -0,250100 -0,250099 -0,250000 -0,240100 -0,210100 -0,160100 -0,090100 -0,000100 0,749900 0,7 -0,490000 -0,489999 -0,489900 -0,480000 -0,450000 -0,400000 -0,330000 -0,240000 0,510000 0,8 -0,640000 -0,639999 -0,639900 -0,630000 -0,600000 -0,550000 -0,480000 -0,390000 0,360000 0,9 -0,810000 -0,809999 -0,809900 -0,800000 -0,770000 -0,720000 -0,650000 -0,560000 0,190000 0,9001 -0,810180 -0,810179 -0,810080 -0,800180 -0,770180 -0,720180 -0,650180 -0,560180 0,189820 0,901 -0,811801 -0,811800 -0,811701 -0,801801 -0,771801 -0,721801 -0,651801 -0,561801 0,188199 0,91 -0,828100 -0,828099 -0,828000 -0,818100 -0,788100 -0,738100 -0,668100 -0,578100 0,171900 1 -1,000000 -0,999999 -0,999900 -0,990000 -0,960000 -0,910000 -0,840000 -0,750000 0,000000

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

79

Tabla 2: Valores numéricos de la solución en la región de la diagonal x=y.

X/Y 0,390000 0,399000 0,399900 0,399990 0,399999 0,400000 0,400001 0,40001 0,4001 0,401 0,41

0 0,152100 0,159201 0,159920 0,159992 0,159999 0,160000 0,160001 0,160008 0,160080 0,160801 0,168100

0,000001 0,152100 0,159201 0,159920 0,159992 0,159999 0,160000 0,160001 0,160008 0,160080 0,160801 0,168100

0,00001 0,152100 0,159201 0,159920 0,159992 0,159999 0,160000 0,160001 0,160008 0,160080 0,160801 0,168100

0,0001 0,152100 0,159201 0,159920 0,159992 0,159999 0,160000 0,160001 0,160008 0,160080 0,160801 0,168100

0,001 0,152099 0,159200 0,159919 0,159991 0,159998 0,159999 0,160000 0,160007 0,160079 0,160800 0,168099

0,01 0,152000 0,159101 0,159820 0,159892 0,159899 0,159900 0,159901 0,159908 0,159980 0,160701 0,168000

0,1 0,142100 0,149201 0,149920 0,149992 0,149999 0,150000 0,150001 0,150008 0,150080 0,150801 0,158100

0,39 0,000000 0,007101 0,007820 0,007892 0,007899 0,007900 0,007901 0,007908 0,007980 0,008701 0,016000 0,399 -0,007101 0,000000 0,000719 0,000791 0,000798 0,000799 0,000800 0,000807 0,000879 0,001600 0,008899

0,3999 -0,007820 -0,000719 0,000000 0,000072 0,000079 0,000080 0,000081 0,000088 0,000160 0,000881 0,008180

0,39999 -0,007892 -0,000791 -0,000072 0,000000 0,000007 0,000008 0,000009 0,000016 0,000088 0,000809 0,008108

0,399999 -0,007899 -0,000798 -0,000079 -0,000007 0,000000 0,000001 0,000002 0,000009 0,000081 0,000802 0,008101

0,4 -0,007900 -0,000799 -0,000080 -0,000008 -0,000001 0,000000 0,000001 0,000008 0,000080 0,000801 0,008100

0,400001 -0,007901 -0,000800 -0,000081 -0,000009 -0,000002 -0,000001 0,000000 0,000007 0,000079 0,000800 0,008099

0,40001 -0,007908 -0,000807 -0,000088 -0,000016 -0,000009 -0,000008 -0,000007 0,000000 0,000072 0,000793 0,008092

0,4001 -0,007980 -0,000879 -0,000160 -0,000088 -0,000081 -0,000080 -0,000079 -0,000072 0,000000 0,000721 0,008020

0,401 -0,008701 -0,001600 -0,000881 -0,000809 -0,000802 -0,000801 -0,000800 -0,000793 -0,000721 0,000000 0,007299

0,41 -0,016000 -0,008899 -0,008180 -0,008108 -0,008101 -0,008100 -0,008099 -0,008092 -0,008020 -0,007299 0,000000

0,411 -0,016821 -0,009720 -0,009001 -0,008929 -0,008922 -0,008921 -0,008920 -0,008913 -0,008841 -0,008120 -0,000821

0,4111 -0,016903 -0,009802 -0,009083 -0,009011 -0,009004 -0,009003 -0,009002 -0,008995 -0,008923 -0,008202 -0,000903

0,41111 -0,016911 -0,009810 -0,009091 -0,009019 -0,009012 -0,009011 -0,009011 -0,009003 -0,008931 -0,008210 -0,000911

0,49 -0,088000 -0,080899 -0,080180 -0,080108 -0,080101 -0,080100 -0,080099 -0,080092 -0,080020 -0,079299 -0,072000

0,499 -0,096901 -0,089800 -0,089081 -0,089009 -0,089002 -0,089001 -0,089000 -0,088993 -0,088921 -0,088200 -0,080901

0,4999 -0,097800 -0,090699 -0,089980 -0,089908 -0,089901 -0,089900 -0,089899 -0,089892 -0,089820 -0,089099 -0,081800

0,49999 -0,097890 -0,090789 -0,090070 -0,089998 -0,089991 -0,089990 -0,089989 -0,089982 -0,089910 -0,089189 -0,081890

0,499999 -0,097899 -0,090798 -0,090079 -0,090007 -0,090000 -0,089999 -0,089998 -0,089991 -0,089919 -0,089198 -0,081899

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

80

3. MÉTODO DE LOS ELEMENTOS FINITOS:

MÉTODO NUMÉRICO ELEGIDO PARA LA RESOLUCIÓN:

Para resolver la ecuación del calor numéricamente se recurre al

método de los elementos finitos.

Este método es aplicable a multitud de sistemas de ecuaciones

diferenciales, de las cuales la Ecuación del Calor no es más que

un caso particular aunque muy ilustrativo.

El método de los elementos finitos (MEF en castellano o FEM en

inglés) es un método numérico general para la aproximación de

soluciones de ecuaciones diferenciales parciales muy utilizado en

diversos problemas de ingeniería y física.

El MEF está pensado para ser usado en computadoras y permite

resolver ecuaciones diferenciales asociadas a un problema físico

sobre geometrías complicadas. El MEF se usa en el diseño y mejora

de productos y aplicaciones industriales, así como en la

simulación de sistemas físicos y biológicos complejos. La

variedad de problemas a los que puede aplicarse ha crecido

enormemente, siendo el requisito básico para su uso el que las

ecuaciones constitutivas y ecuaciones de evolución temporal del

problema a considerar sean conocidas de antemano.

En las figuras 1 y 2 se ilustra un ejemplo de resolución mediante

el uso de este método.

Fig. 1 Solución de MEF en 2D para una configuración de

un magnetostato, (las líneas muestran la dirección de la

densidad de flujo calculada, y el color, su magnitud).

Fig. 2 La malla 2D para la imagen superior (la malla es más

densa alrededor del nuestro objetivo, aquellas zonas de mayor interés, o de mayor complejidad en el cálculo).

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

81

El MEF permite obtener una solución numérica aproximada sobre un cuerpo,

estructura o dominio (medio continuo) sobre el que están definidas

ciertas ecuaciones diferenciales en forma débil o integral que

caracterizan el comportamiento físico del problema dividiéndolo en un

número elevado de subdominios no-intersectantes entre sí denominados

«elementos finitos». El conjunto de elementos finitos forma una

partición del dominio también denominada discretización. Dentro de cada

elemento se distinguen una serie de puntos representativos llamados

«nodos». Dos nodos son adyacentes si pertenecen al mismo elemento

finito; además, un nodo sobre la frontera de un elemento finito puede

pertenecer a varios elementos. El conjunto de nodos considerando sus

relaciones de adyacencia se llama «malla».

Los cálculos se realizan sobre una malla de puntos (llamados nodos), que

sirven a su vez de base para la discretización del dominio en elementos

finitos. La generación de la malla se realiza usualmente con programas

especiales llamados generadores de mallas, en una etapa previa a los

cálculos que se denomina pre-proceso. De acuerdo con estas relaciones de

adyacencia o conectividad se relaciona el valor de un conjunto de

variables incógnitas definidas en cada nodo y denominadas grados de

libertad. El conjunto de relaciones entre el valor de una determinada

variable entre los nodos se puede escribir en forma de un sistema de

ecuaciones lineales (o linealizadas). La matriz de dicho sistema de

ecuaciones se llama matriz de rigidez del sistema. El número de

ecuaciones de dicho sistema es proporcional al número de nodos.

Habitualmente, para el caso de un problema de Mecánica de los Sólidos

Deformables o más generalmente un problema de Mecánica de Medios

Continuos, el método de los elementos finitos se programa

computacionalmente, para calcular el campo de desplazamientos y,

posteriormente, a través de relaciones cinemáticas y constitutivas, las

deformaciones y tensiones respectivamente. El método de los elementos

finitos es muy usado debido a su generalidad y a la facilidad de

introducir dominios de cálculo complejos (en dos o tres dimensiones).

Además el método es fácilmente adaptable a problemas de transmisión de

calor, de mecánica de fluidos para calcular campos de velocidades y

presiones (mecánica de fluidos computacional, CFD) o de campo

electromagnético. Dada la imposibilidad práctica de encontrar la

solución analítica de estos problemas, con frecuencia en la práctica

ingenieril los métodos numéricos y, en particular, los elementos

finitos, se convierten en la única alternativa práctica de cálculo.

Una importante propiedad del método es la convergencia; si se consideran

particiones de elementos finitos sucesivamente más finas, la solución

numérica calculada converge rápidamente hacia la solución exacta del

sistema de ecuaciones.

PROGRAMACIÓN EN EL ENTORNO CUDA EN APLICACIONES DE MECÁNICA COMPUTACIONAL

CÓDIGO DEL PROYECTO

82

El desarrollo de un algoritmo de elementos finitos para resolver un

problema definido mediante ecuaciones diferenciales y condiciones de

contorno requiere en general cuatro etapas:

El problema debe reformularse en forma variacional.

El dominio de variables independientes (usualmente un dominio

espacial para problemas dependientes del tiempo) debe dividirse

mediante una partición en subdominios, llamados elementos finitos.

Asociada a la partición anterior se construye un espacio vectorial

de dimensión finita, llamado espacio de elementos finitos. Siendo

la solución numérica aproximada obtenida por elementos finitos una

combinación lineal en dicho espacio vectorial.

Se obtiene la proyección del problema variacional original sobre

el espacio de elementos finitos obtenido de la partición. Esto da

lugar a un sistema con un número de ecuaciones finito, aunque en

general con un número elevado de ecuaciones incógnitas. El número

de incógnitas será igual a la dimensión del espacio vectorial de

elementos finitos obtenido y, en general, cuanto mayor sea dicha

dimensión tanto mejor será la aproximación numérica obtenida.

El último paso es el cálculo numérico de la solución del sistema

de ecuaciones.

Los pasos anteriores permiten convertir un problema de cálculo

diferencial en un problema de álgebra lineal. Dicho problema en

general se plantea sobre un espacio vectorial de dimensión no-

finita, pero que puede resolverse aproximadamente encontrando una

proyección sobre un espacio subespacio de dimensión finita, y por

tanto con un número finito de ecuaciones (aunque en general el

número de ecuaciones será elevado típicamente de miles o incluso

centenares de miles). La discretización en elementos finitos

ayuda a construir un algoritmo de proyección sencillo, logrando

además que la solución por el método de elementos finitos sea

generalmente exacta en un conjunto finito de puntos. Estos puntos

coinciden usualmente con los vértices de los elementos finitos o

puntos destacados de los mismos. Para la resolución concreta del

enorme sistema de ecuaciones algebraicas en general pueden usarse

los métodos convencionales del álgebra lineal en espacios de

dimensión finita.