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.