Source: http://docplayer.es/1501906-Metodo-de-descomposicion-de-dominio-de-subestructuracion.html
Timestamp: 2016-10-27 03:41:45+00:00

Document:
⭐Método de Descomposición de Dominio de Subestructuración
Método de Descomposición de Dominio de Subestructuración
Download "Método de Descomposición de Dominio de Subestructuración"
Patricia Pérez Duarte
1 Método de Descomposición de Dominio de Subestructuración Antonio Carrillo Ledesma Febrero 2012 Índice 1. Introducción MétodosdeDescomposicióndeDominio ObjetivosdelTrabajo ObjetivosGenerales ObjetivosParticulares InfraestructuraUsada Análisis Funcional y Problemas Variacionales OperadorLinealElíptico EspaciosdeSobolev TrazasdeunaFunciónenH m () EspaciosH m 0 () ProblemasVariacionalesconValorenlaFrontera Solución de Grandes Sistemas de Ecuaciones MétodosDirectos DescomposiciónLU DescomposicióndeCholesky MétodosIterativos MétododeGradienteConjugado MétodoResidualMínimoGeneralizado Estructura óptima de las Matrices en su Implementación Computacional MatricesBandadas MatricesDispersas MultiplicaciónMatriz-Vector Métodos de Solución Aproximada para EDP MétodoGalerkin ElMétododeResiduosPesados MétododeElementosFinitos2 5. Método de Elementos Finitos Triangulación InterpolaciónparaelMétododeElementosFinitos Método de Elemento Finito Usando Discretización de Rectángulos Método de Elemento Finito Usando Discretización de Triángulos Métodos de Descomposición de Dominio (DDM) MétododeSchwarz Método de Descomposición de Dominio de Subestructuración El Cómputo en Paralelo ArquitecturasdeSoftwareyHardware ClasificacióndeFlynn CategoríasdeComputadorasParalelas MétricasdeDesempeño CómputoParaleloparaSistemasContinuos Implementación Computacional Secuencial y Paralela de DDM ElOperadordeLaplaceylaEcuacióndePoisson MétododelElementoFinitoSecuencial MétododeSubestructuraciónSecuencial MétododeSubestructuraciónenParalelo Método de Subestructuración en Paralelo Precondicionado Análisis de Rendimiento y Conclusiones AnálisisdeComunicaciones Afectación del Rendimiento al Aumentar el Número de SubdominiosenlaDescomposición DescomposiciónÓptimaparaunEquipoParaleloDado ConsideracionesparaAumentarelRendimiento ConclusionesGenerales Apéndice A NocionesdeAlgebraLineal σ-AlgebrayEspaciosMedibles EspaciosL p Distribuciones Bibliografía 119 23 1. Introducción La necesidad de entender su entorno y anticiparse a los acontecimientos tiene raíces muy profundas en el ser humano. Desde la prehistoria, el hombre trató de predecir a la naturaleza, pues de ella dependía su supervivencia, para lo cual inicialmente nuestros antepasados utilizaron a la brujería, así como el pensamiento mágico y el religioso. Sin embargo, el medio más efectivo para predecir el comportamiento de la naturaleza es el método científico(o su antecesor el método empírico)yesporesoqueesteanhelohumanoancestral,atravésdelahistoria, hasidomotordelaciencia. Lamaduraciónyelprogresodelaprediccióncientíficaes,sinduda,elresultado del avance general de la ciencia, pero además ha habido elementos catalizadores esenciales sin los cuales esto no hubiera sido posible. La predicción científica, además de ser científica, es matemática y computacional. En la actualidad, cuando deseamos predecir el comportamiento de un sistema, los conocimientos científicos y tecnológicos se integran en modelos matemáticos los cuales se convierten en programas de cómputo que son ejecutados por las computadoras tanto secuenciales como paralelas(entenderemos por una arquitectura paralela a un conjunto de procesadores interconectados capaces de cooperar en la solución de un problema). Aunque el sólo hecho de poseer la capacidad de predicción científica nos llena de satisfacción y orgullo, sin embargo, aún más trascendente es el hecho de que ella también es la base de una gran parte del extraordinario progreso material que la humanidad ha experimentado en épocas recientes. En efecto, nuestra facultad para predecir es una herramienta muy poderosa de la ingeniería, de la tecnología y de la ciencia misma, la cual, entre otras muchas cosas, nos ha permitido ampliar la disponibilidad de los recursos naturales y utilizarlos con mayor eficiencia. Para resolver, por ejemplo, algún problema del área de ciencias e ingeniería (elmovimientodeunfluidolibreodeunfluidoenunmedioporoso)loprimero que debemos de hacer es formular un modelo matemático del problema a tratar. Estoselograatravésdelosfundamentosdelafísicamacroscópicaloscualesson proporcionados por la teoría de los medios continuos [?]. Con base en ella se introduce una formulación clara, general y sencilla de los modelos matemáticos de los sistemas continuos. Laformulaciónestansencillaytangeneral,quelosmodelosbásicosdesistemas tan complicados y diversos como la atmósfera, los océanos, los yacimientos petroleros, o los geotérmicos, se derivan por medio de la aplicación repetida de una sola ecuación diferencial: la ecuación diferencial de balance [?]. Dichaformulacióntambiénesmuyclara,puesenelmodelogeneralnohay ninguna ambigüedad; en particular, todas las variables y parámetros que intervienen en él, están definidos de manera unívoca. En realidad, este modelo general de los sistemas continuos constituye una realización extraordinaria de los paradigmas del pensamiento matemático. El descubrimiento del hecho de que los modelos matemáticos de los sistemas continuos, independientemente de su naturaleza y propiedades intrínsecas, pueden formularse por medio de balances, 34 cuya idea básica no difiere mucho de los balances de la contabilidad financiera, fue el resultado de un largo proceso de perfeccionamiento en el que concurrieron una multitud de mentes brillantes. Los modelos matemáticos de los sistemas continuos son ecuaciones diferenciales, las cuales son parciales (con valores iniciales y condiciones de frontera) paracasitodoslossistemasdemayorinterésenlacienciaylaingeniería,osistemas de tales ecuaciones. Salvo para los problemas más sencillos, no es posible obtener por métodos analíticos las soluciones de tales ecuaciones, que son las que permiten predecir el comportamiento de los sistemas. La capacidad para formular los modelos matemáticos de sistemas complicados y de gran diversidad, es sin duda una contribución fundamental para el avance de la ciencia y sus aplicaciones, tal contribución quedaría incompleta y, debido a ello, sería poco fecunda, si no se hubiera desarrollado simultáneamente su complemento esencial: los métodos matemáticos y la computación electrónica. En cambio, la diversidad y complejidad de problemas que pueden ser tratados con métodos numéricos y computacionales es impresionante. Los modelos de los sistemas continuos-es decir, sistemas físicos macroscópicos tales como los yacimientos petroleros, la atmósfera, los campos electromagnéticos, los océanos, los metalúrgicos, el aparato circulatorio de los seres humanos, la corteza terrestre, lo suelos y las cimentaciones, muchos sistemas ambientales, y muchos otros cuya enumeración ocuparía un espacio enormecontienen un gran número de grados libertad. Por ello, la solución numérica por los esquemas tradicionales tipo diferencias finitas y elemento finito generan una discretización del problema, la cual es usada para generar sistemas de ecuaciones algebraicos[?]. Estos sistemas algebraicos en general son de gran tamaño para problemas reales, al ser estos algoritmos secuenciales su implantación suele hacerse en equipos secuenciales y por ello no es posible resolver muchos problemas que involucren el uso de una gran cantidad de memoria. Actualmente para tratar de subsanar la limitante de procesar sólo en equipos secuenciales, se usan equipos paralelos para soportar algoritmos secuenciales mediante directivas de compilación, haciendo ineficiente su implantación en dichos equipos. La computación en paralelo es una técnica que nos permite distribuir una gran carga computacional entre muchos procesadores. Y es bien sabido que una de las mayores dificultades del procesamiento en paralelo es la coordinación de las actividades de los diferentes procesadores y el intercambio de información entre los mismos[12]. Al usar métodos de descomposición de dominio conjuntamente con el cómputo en paralelo(supercomputadoras, clusters o grids) nos permite atacar una gran variedad de problemas que sin estas técnicas sería imposible hacerlo de manera eficiente Métodos de Descomposición de Dominio Los métodos de descomposición de dominio introducen desde la formulación matemática del problema una separación natural de las tareas a realizar por el 45 método y simplifican considerablemente la transmisión de información entre los subdominios[?]. En ellos, los sistemas físicos representados por su modelo son descompuestos en varios subdominios. Los métodos de descomposición de dominio permiten tratar los problemas de tamaño considerable, empleando algoritmos paralelos en computadoras secuenciales y/o paralelas. Esto es posible ya que cualquier método de descomposición dedominiosebasaenlasuposicióndequedadoundominiocomputacional, estesepuedeparticionarensubdominios i, i=1,2,...,eentreloscualespuede o no existir traslape[1] y[?]. Entonces, el problema es reformulado en términos de cada subdominio (empleando algún método del discretización) obteniendo una familia de subproblemas de tamaño reducido independientes en principio entre si, y que están acoplados a través de la solución en la interfaz de los subdominios que es desconocida. De esta manera, podemos clasificar de forma burda a los métodos de descomposición de dominio, como aquellos en que: existe traslape entre los subdominiosyenlosquenoexistetraslape.alaprimeraclaseperteneceelmétodo deschwarz(enelcualeltamañodeltraslapeesimportanteenlaconvergencia del método) y a los de la segunda clase pertenecen los métodos del tipo subestructuración(en el cual los subdominios sólo tienen en común los nodos de la frontera interior). Estos métodos son un paradigma natural usado por la comunidad de modeladores. Los sistemas físicos son descompuestos en dos o más subdominios contiguos basados en consideraciones fenomenológicas. Estas descomposiciones basadas en dominios físicos son reflejadas en la ingeniería de software del código correspondiente. Así, mediante los métodos de descomposición de dominio, la programación orientada a objetos (que nos permite dividir en niveles la semántica de los sistemas complejos tratando así con las partes, más manejables que el todo, permitiendo su extensión y un mantenimiento más sencillo) y esquemas de paralelización que usan el paso de mensajes, es posible construir aplicaciones que coadyuven a la solución de problemas concomitantes en ciencia e ingeniería. Esta metodología permite utilizar todas las capacidades del cómputo en paralelo(grids de decenas de clusters, cada uno con miles de procesadores interconectados por red con un creciente poder de cómputo medible en Peta Flops), asícomoelusodeunaampliamemoria(yaseadistribuiday/ocompartidadel orden de Tera Bytes), permitiendo atacar una gran variedad de problemas que sin estas técnicas es imposible hacerlo de manera flexible y eficiente. Pero hay que notar que existe una amplia gama de problemas que nos interesan resolver que superan las capacidades de cómputo actuales, ya sea por el tiempo requerido para su solución, por el consumo excesivo de memoria o ambos. Lalistadelosmétodosdedescomposicióndedominioyeltipodeproblemas que pueden ser atacados por estos, es grande y está en constante evolución[?], ya que se trata de encontrar un equilibrio entre la complejidad del método (aunadaalapropiacomplejidaddelmodelo),laeficienciaenelconsumodelos recursos computacionales y la precisión esperada en la solución encontrada por 56 los diversos métodos y las arquitecturas paralelas en la que se implante Objetivos del Trabajo Unodelosgrandesretosdeláreadecómputocientíficoespoderanalizara priori una serie de consideraciones dictadas por factores externos al problema de interés que repercuten directamente en la forma de solucionar el problema, estas consideraciones influirán de manera decisiva en la implementación computacional de la solución numérica. Algunas de estas consideraciones son: Número de Procesadores Disponibles TamañoyTipodeParticióndelDominio Tiempo de Ejecución Predeterminado Siendo común que ellas interactúan entre si, de forma tal que normalmente el encargado de la implementación computacional de la solución numérica tiene además de las complicaciones técnicas propias de la solución, el conciliarlas con dichas consideraciones. Esto deja al implementador de la solución numérica con pocos grados de libertad para hacer de la aplicación computacional una herramienta eficiente y flexible que cumpla con los lineamientos establecidos a priori y permita también que esta sea adaptable a futuros cambios de especificaciones -algo común en ciencia e ingeniería-. Paratratardetenerunaideaclaradecómoafectanalaaplicacióncomputacional estos factores, a continuación detallamos algunas consideraciones acerca decadaunodeellos. Número de Procesadores Disponibles El número de procesadores disponibles es una barrera en constante evolución pero las principales limitantes son económicas y tecnológicas, es común ahora contar con miles de procesadores interactuando de manera conjunta en grids, pero todo ese poder de cómputo creciente es aun insuficiente para las necesidades del cómputo científico. Por otro lado la gran mayoría de los proyectos científicos y tecnológicos no cuentan con recursos computacionales grandes, ya que este es un recurso costoso y por ende muy limitado y en la mayoría de los casos no es el adecuado a las necesidades propias del problema. Así, para poder atacar el problema de forma eficiente es necesario adaptarse al equipo de cómputo disponible y tratar de conocer de antemano los factores que mermaran el rendimiento de la implementación computacional para buscar alternativas que mejoren la eficiencia de la implementación. Tamaño y Tipo de Partición del Dominio Normalmente cuando se plantea un problema de sistemas continuos la determinación del dominio y tipo de malla a usar en la solución del mismo esta sujeta a restricciones fenomenológicas, por ello la malla deberá de adaptarse de la mejor forma posible para tratar 67 decapturarlosrasgosesencialesdelfenómenoestudiado. Peroalahoradela implementación computacional es común hacer adecuaciones a esta, para que puedasersoportadaporelequipodecómputoysuejecuciónseaenuntiempo razonable. Esto puede ser un problema sobre todos al ser implementada la solución usando cómputo paralelo, ya que comúnmente una elección de una malla no homogénea ocasionará problemas de mal balanceo de carga de trabajo entre los procesadores utilizados en la implementación computacional. Tiempo de Ejecución Predeterminado Otro factor determinante es el tiempo de solución esperado de la implementación computacional del problema, estoescríticoenproblemasdecontrolentiemporeal,comunesenlacienciae ingeniería actual. Por ello en estos casos es permisible el aumento en el número y capacidades del equipo de cómputo necesario en la implementación con el fin de lograr un tiempo de ejecución por debajo del máximo permisible dictado por las especificaciones propias del problema. Todos estos factores influirán en la versión final implementada en la solución computacional de un problema particular, debiendo ser todas ellas sopesadas antesdetomarunadecisiónencuantoalasespecificacionesqueelprogramade cómputo deberá satisfacer. Así, el objetivo de este trabajo es desarrollar una metodología en la cual se implemente de forma paralela el método numérico de descomposición de dominio de subestructuración precondicionado, explicando detalladamente los fundamentos matemáticos de la modelación de fenómenos de sistemas continuos por medio de ecuaciones diferenciales parciales, los fundamentos del método de descomposición de dominio al aplicarse a ecuaciones diferenciales parciales elípticas y las ventajas sobre otros métodos de solución numérica. Mostraremos como construir el modelo computacional, esto no sólo nos ayudará a demostrar que es factible la construcción del propio modelo computacional a partir del modelo matemático y numérico para la solución de problemas reales. Además, mostrará los alcances y limitaciones en el consumo de recursos computacionales y nos permitirá la evaluación de algunas de las variantes de los métodos numéricos con los que es posible implementar el modelo computacional y haremos algo de análisis de rendimiento sin llegar a ser exhaustivo esté. También se muestra el diseño, análisis y programación de la aplicación computacional del método de descomposición de dominio en forma serial como paralela, basada en el paradigma de programación orientado a objetos y como este paradigma ofrece una forma robusta y flexible para la representación de entidades abstractas y que permiten una modelación más versátil así como una forma más eficiente para la organización y reutilización de código Objetivos Generales Mostrarlasbasesdeunametodologíaqueseutilizaparaaplicarelcómputo en paralelo a la modelación matemática y computacional de sistemas continuos, de forma flexible, escalable y eficiente. 78 Mostrar los alcances y limitaciones de la metodología usando como herramientas de evaluación a los métodos de elemento finito secuencial, método de subestructuración secuencial y método de subestructuración paralelo. Mostrar los diversos esquemas de optimización de los recursos computacionales aplicables a un problema especifico. Para facilitar la comprensión de las ideas básicas, se ha tomado la ecuación de Poisson. Es un ejemplo sencillo, pero gobierna los modelos de muchos sistemas delaingenieríaydelaciencia,entreelloselflujodeaguasubterráneaatravés de un acuífero isotrópico, homogéneo bajo condiciones de equilibrio y es muy usadaenmúltiplesramasde lafísica.porejemplo, gobiernalaecuacióndela conducción de calor en un sólido bajo condiciones de equilibrio Objetivos Particulares Como objetivos particulares de este trabajo tenemos: Mostrar cómo aplicar la metodología para manejar problemas de gran tamaño(descomposición de malla fina). Mostrar cómo descomponer un dominio en un conjunto de subdominios quedenunaparticiónenlaqueeltiempodecálculoseamínimoparauna configuración de hardware dada. Mostrar cuales son las posibles optimizaciones aplicables a una configuración de hardware dada. Mostrarqueesposibletrabajarproblemasconunamallamuyfinaenun equipo paralelo pequeño. Para poder cumplir con los objetivos, primeramente empezaremos describiendo las bases de los sistemas continuos y sus modelos, para conocer algunas propiedades de las ecuaciones diferenciales parciales que son generadas en la elaboración del modelo matemático. Después veremos la forma de pasar del modelomatemáticoalmodelonumérico(elcualnoesmuycercanoalmodelocomputacional) introduciendo los fundamentos del método Galerkin y la derivación a partir de este del método de elemento finito, para posteriormente describir dos métodos de descomposición de dominio(schwarz y subestructuración). Mostraremos como construir el modelo computacional que dependerá fuertemente de la arquitectura de cómputo disponible y del lenguaje de programación usado. Esto nos ayudará a demostrar que es factible la construcción del propio modelo computacional a partir del modelo matemático y numérico para la solución de problemas reales. Además conoceremos los alcances y limitaciones en el consumo de recursos computacionales y nos permitirá la evaluación de algunas de las variantes de los 89 métodos numéricos con los que es posible implementar el modelo computacional y haremos algo de análisis de rendimiento sin llegar a ser exhaustivo este. Finalmente se exploran los alcances y limitaciones de cada uno de los métodos implementados(elemento finito secuencial, descomposición de dominio secuencial y paralelo) y se muestra como es posible optimizar los recursos computacionales con que se cuente Infraestructura Usada El modelo computacional generado, está contenido en un programa de cómputo bajo el paradigma de orientación a objetos, programado en el lenguaje C++ en su forma secuencial y en su forma paralela en C++ en conjunto con la interfaz de paso de mensajes (MPI) bajo el esquema de maestroesclavo. Una versión de estos programas está disponible en la página Web bajo el rubro FEM y DDM. Para desarrollar estos códigos, se realizo una jerarquía de clases para cada uno de los distintos componentes del sistema de elemento finito como de descomposición de dominio(rehusando todo el código de elemento finito en este), permitiendo usarlos tanto en forma secuencial como paralela. Las pruebas de rendimiento de los distintos programas se realizaron en equipos secuenciales y paralelos (clusters) que están montados en el Instituto degeofísicadelaunam,enlaspruebasdeanálisisderendimientoseusaron para la parte secuencial el equipo: Computadora Pentium IV HT a 2.8 GHtz con 1 GB de RAM corriendo bajo el sistema operativo Linux Debian Stable con el compilador g++ de GNU. Para la parte paralela se usaron los siguientes equipos: Cluster heterogéneo con el nodo maestro Pentium IV HT a 3.4 GHtzcon1GBdeRAMy7nodosesclavosPentiumIVHTa2.8 GHtz con 0.5 GB de RAM por nodo, unidos mediante una red E- thernet de 100 Mb, corriendo bajo el sistema operativo Linux Debian StableconelcompiladormpiCCdeMPIdeGNU.AcargodelDr. Ismael Herrera Revilla del departamento de Recursos Naturales. Hay que notar que, el paradigma de programación orientada a objetos sacrifica algo de eficiencia computacional por requerir mayor manejo de recursos computacionales al momento de la ejecución. Pero en contraste, permite mayor flexibilidad a la hora adaptar los códigos a nuevas especificaciones. Adicionalmente, disminuye notoriamente el tiempo invertido en el mantenimiento y búsqueda de errores dentro del código. Esto tiene especial interés cuando se piensa en la cantidad de meses invertidos en la programación comparado con los segundos consumidos en la ejecución del mismo. 910 2. Análisis Funcional y Problemas Variacionales En esta sección se detallan los conceptos básicos de análisis funcional y problemas variacionales con énfasis en problemas elípticos de orden par 2m, para comenzar detallaremos lo que entendemos por un operador diferencial parcial elípticodeordenpar2mennvariables,paradespuésdefiniralosespaciosde Sobolev para poder tratar problemas variacionales con valor en la frontera. En donde, restringiéndonos a problemas elípticos, contestaremos una cuestión centralenlateoríadeproblemaselípticosconvaloresenlafrontera,yestáse relaciona con las condiciones bajo las cuales uno puede esperar que el problema tenga solución y esta es única, así como conocer la regularidad de la solución, paramayorreferenciadeestosresultadosver[?],[?]y[?] Operador Lineal Elíptico Definición1 Entenderemosporundominioalconjunto R n queseaabiertoyconexo. Para poder expresar de forma compacta derivadas parciales de orden m o menor, usaremos la definición siguiente. Definición2 Sea Z n + el conjunto de todas las n-duplas de enteros no negativos, un miembro de Z n + se denota usualmente por α ó β (por ejemplo α = (α 1,α 2,...,α n ).Denotaremospor α lasuma α =α 1 +α α n ypord α u la derivada parcial D α α u u= (1) x α1 1 xα xαn n así,si α =m,entoncesd α udenotalam-ésimaderivadaparcialdeu. SeaLunoperadordiferencialparcialdeordenpar2mennvariablesyde la forma Lu= ( 1) α D α( a αβ (x)d β u ), x R n (2) α, β m donde es un dominio en R n. Los coeficientes a αβ son funciones suaves real valuadas de x. El operador L es asumido que aparece dentro de una ecuación diferencial parcial de la forma Lu=f, (3) dondef pertenecealrangodeloperadorl. La clasificación del operador L depende sólo de los coeficientes de la derivada másalta,estoes,deladerivadadeorden2m,yalostérminosinvolucradosen esaderivadasonllamadoslaparteprincipaldeloperadorldenotadoporl 0 y paraeloperador(2)esdelaforma L 0 = a αβ D α+β u. (4) α, β m 1011 Teorema3 Seaξ unvectorenr n,yseaξ α =ξ α1 n...ξαn n,α Z+ n.entonces i)leselípticoenx o si α, β =m a αβ ( xo ) ξ α+β 0 ξ 0; (5) ii)leselípticosieselípticoentodoslospuntosde; iii)lesfuertementeelípticosiexisteunnúmeroµ>0talque ( ) a αβ xo ξ α+β µ ξ 2m α, β =m (6) satisfaciéndoseentodopuntox o,yparatodoξ R n.aquí ξ = ( ξ ξ 2 n)1 2. ParaelcasoenelcualLesunoperadorde2doorden(m=1),lanotación se simplifica, tomando la forma Lu= n i,j=1 x i ( a ij (x) u ) x j + n j=1 a j u x j +a 0 u=f (7) en. Paracoeficientesadecuadosa ij,a j ya 0 lacondiciónparaconocersieloperador es elíptico, es examinado por la condición n a ij (x 0 )ξ i ξ j 0 ξ 0 (8) i,j=1 y para conocer si el operador es fuertemente elíptico, es examinado por la condición n a ij (x 0 )ξ i ξ j >µ ξ 2. (9) i,j= Espacios de Sobolev En esta subsección detallaremos algunos resultados de los espacios de Sobolev sobre el conjunto de números reales, en estos espacios son sobre los cuales trabajaremos tanto para plantear el problema elíptico como para encontrar la solución alproblema.primeramentedefiniremosloqueentendemosporunespaciol 2. Definición4 Una función medible u(x) definida sobre R n se dice que pertenecealespaciol 2 ()si u(x) 2 dx< (10) es decir, es integrable. 1112 Ladefinicióndelosespaciosmedibles,espaciosL p,distribucionesyderivadas de distribuciones están dados en el apéndice, estos resultados son la base para poder definir a los espacios de Sobolev. Definición5 ElespaciodeSobolevdeordenm,denotadoporH m (),esdefinido H m ()= { u D α u L 2 () αtalque α m }. (11) Elproductoescalar, dedoselementosuyv H m ()estadadopor u,v H m= (D α u)(d α v)dxparau,v H m (). (12) α m Nota:EscomúnqueelespacioL 2 ()seadenotadoporh 0 (). Un espacio completo con producto interior es llamado un espacio de Hilbert, un espacio normado y completo es llamado espacio de Banach. Y como todo producto interior define una norma, entonces todo espacio de Hilbert es un espacio de Banach. Definición6 Lanorma H m inducidaapartirdelproductointerior, H m queda definida por u 2 H m= u,u H m= (D α u) 2 dx. (13) α m Ahora, connorma H m,el espacioh m ()es unespaciodehilbert, esto queda plasmado en el siguiente resultado. Teorema7 ElespacioH m ()conlanorma H m esunespaciodehilbert. YaquealgunasdelaspropiedadesdelosespaciosdeSobolevsólosonvalidas cuando la frontera del dominio es suficientemente suave. Para describir al conjunto donde los espacios de Sobolev están definidos, es común pedirle algunas propiedades y así definimos lo siguiente. Definición8 Una función f definida sobre un conjunto Γ R n es llamada LipschitzcontinuasiexisteunaconstanteL>0talque f(x) f(y) L x y x,y Γ. (14) Notemos que una función Lipschitz continua es uniformemente continua. Sea R n (n 2) un dominio con frontera, sea x 0 y construyamos la bola abierta con centro en x 0 y radio ε, i.e. B(x 0,ε), entonces definiremoselsistemacoordenado(ξ 1,...,ξ n )talqueelsegmento B(x 0,ε) pueda expresarse como una función entonces definimos. ξ n =f(ξ 1,...,ξ n 1 ) (15) 1213 Definición9 Lafrontera deldominioesllamadadelipschitzsif definidacomoenlaec.(15)esunafunciónlipschitzcontinua. El siguiente teorema resume las propiedades más importantes de los espacios desobolevh m (). Teorema10 Sea H m () el espacio de Sobolev de orden m y sea R n un dominio acotado con frontera Lipschitz. Entonces i)h r () H m ()sir m ii)h m ()esunespaciodehilbertconrespectoalanorma H m iii)h m ()eslacerraduraconrespectoalanorma H mdelespacioc (). De la parte iii) del teorema anterior, se puede hacer una importante interpretación: Para toda u H m () es siempre posible encontrar una función infinitamente diferenciable f, tal que este arbitrariamente cerca de u en el sentido que u f H m<ε (16) paraalgúnε>0dado. Cuandom=0,sededucelapropiedadH 0 ()=L 2 ()apartirdelteorema anterior. Corolario 11 El espacio L 2 () es la cerradura, con respecto a la norma L 2, delespacioc (). Otrapropiedad,setienealconsideraracualquiermiembrodeu H m (), estepuedeseridentificadoconunafunciónenc m (),despuésdequeposiblemente sean cambiados algunos valores sobre un conjunto de medida cero, esto queda plasmado en los dos siguientes resultados. Teorema12 SeanX yy dosespaciosdebanach,conx Y.Seai:X Y tal que i(u)=u. Si el espacio X tiene definida la norma X y el espacio Y tiene definida la norma Y, decimos que X está inmersa continuamente en Y si i(u) Y = u Y K u X (17) paraalgunaconstantek>0. Teorema 13 (Inmersión de Sobolev) Sea R n undominioacotadoconfrontera delipschitz.si(m k)> n/2,entoncestodafunciónenh m ()perteneceac k (),esdecir,hayunmiembroqueperteneceac k ().Además,lainmersión es continua. H m () C k () (18) 1314 Trazasde unafunciónenh m (). Una parte fundamental en los problemas con valores en la frontera definidos sobreeldominio,esdefinirdeformaúnicalosvaloresquetomarálafunción sobre la frontera, en este apartado veremos bajo que condiciones es posible tener definidos de forma única los valores en la frontera tal que podamos definirunoperadortr( )continuoqueactuéentalquetr(u)=u. Elsiguientelemanosdicequeeloperadortr( )esunoperadorlinealcontinuo dec 1( ) ac( ),conrespectoalasnormas H1 () y L 2 ( ). Lema 14 Seaundominioconfrontera delipschitz.laestimación tr(u) L2 ( ) C u H 1 () (19) sesatisfaceparatodafunciónu C 1( ),paraalgunaconstantec>0. Ahora,paraelcasotr( ):H 1 () L 2 ( ),setieneelsiguienteteorema. Teorema15 Sea un dominio acotado en R n con frontera de Lipschitz. Entonces: i)existeunúnicooperadorlinealacotadotr( ):H 1 () L 2 ( ),talque tr(u) L 2 ( ) C u H 1 (), (20) conlapropiedadquesiu C 1( ),entoncestr(u)=u. ii)elrangodetr( )esdensoenl 2 ( ). ElargumentoanteriorpuedesergeneralizadoparalosespaciosH m (),de hecho,cuandom>1,entoncesparatodau H m ()tenemosque D α u H 1 () para α m 1, (21) porelteoremaanterior,elvalorded α usobrelafronteraestábiendefinidoy perteneceal 2 (),esdecir tr(d α u) L 2 (), α m 1. (22) Además,siuesm-vecescontinuamentediferenciable,entoncesD α uesalmenos continuamentediferenciablepara α m 1y Espacios H m 0 (). tr(d α u)=(d α u). (23) LosespacioH m 0 ()surgencomúnmentealtrabajarconproblemasconvalor en la frontera y serán aquellos espacios que se nulifiquen en la frontera del dominio, es decir 1415 Definición16 DefinimosalosespaciosH0 m()comolacerradura,enlanor- madesobolev H m,delespacioc0 m ()defuncionesconderivadascontinuas del orden menor quem,todas las cuales tienen soporte compacto en, es decir H0 m() es formado al tomar la unión de Cm 0 () y de todos los límites de sucesionesdecauchyenc0 m ()quenopertenecenac0 m (). Las propiedades básicas de estos espacios están contenidas en el siguiente resultado. Teorema17 SeaundominioacotadoenR n confrontera suficientemente suaveyseah m 0 ()lacerraduradec 0 ()enlanorma H m,entonces a)h m 0 ()eslacerraduradec 0 ()enlanorma H m; b)h m 0 () Hm (); c)siu H m ()perteneceah m 0 (),entonces D α u=0, sobre, α m 1. (24) Teorema 18 (Desigualdad de Poincaré-Friedrichs) SeaundominioacotadoenR n.entoncesexisteunaconstantec>0tal que u 2 dx C u 2 dx (25) paratodau H 1 0(). Introduciendo ahora una familia de semi-normas sobre H m () (una seminorma satisfacecasitodoslosaxiomasdeunanormaexceptoeldepositivo definido), de la siguiente forma: Definición19 Lasemi-norma m sobreh m (),sedefinecomo u 2 m = D α u 2 dx. (26) α =m Estaes unasemi-norma,ya que u m =0 implicaqued α u=0 para α =m, locualnoimplicaqueu=0. La relevancia de esta semi-norma está al aplicar la desigualdad de Poincaré- Friedrichs ya que es posible demostrar que 1 es de hecho una norma sobre H 1 0(). Corolario 20 Lasemi-norma 1 esunanormasobreh 1 0 (),equivalenteala normaestándar H 1. EsposibleextenderelteoremaanteriorysucorolarioalosespaciosH m 0 () paracualquierm 1,delasiguienteforma: 1516 Teorema21 SeaundominioacotadoenR n.entoncesexisteunaconstante C>0talque u 2 L 2 C u 2 m (27) paratodau H m 0 (),además, m esunanormasobreh m 0 ()equivalentea lanormaestándar H m. Definición22 SeaundominioacotadoenR n.definimosporh m ()alespaciodetodaslasfuncionaleslinealesacotadassobreh m 0 (),esdecir,h m () seráelespaciodualdelespacioh m 0 (). Teorema23 q será una distribución de H m () si y sólo si q puede ser expresada en la forma q= D α q α (28) dondeq α sonfuncionalesenl 2 (). α <m 2.3. Problemas Variacionales con Valor en la Frontera Restringiéndonos ahora en problemas elípticos de orden 2 (problemas de orden mayor pueden ser tratados de forma similar), reescribiremos este en su forma variacional. La formulación variacional es más débil que la formulación convencional ya que esta demanda menor suavidad de la solución u, sin embargo cualquier problema variacional con valores en la frontera corresponde a un problema con valor en la frontera y viceversa. Además, la formulación variacional facilita el tratamiento de los problemas al usar métodos numéricos de ecuaciones diferenciales parciales, en esta sección veremos algunos resultados clave como es la existencia y unicidad de la solución deestetipodeproblemas,paramayoresdetalles,ver[?]y[?]. SieloperadorLestádefinidopor Lu= a u+cu (29) conaunamatrizpositivadefinida,simétricayc 0,elproblemaquedaescrito como a u+cu = f en (30) u = g en. Simultiplicamosalaecuación a u+cu=f porv V =H 1 0(), obtenemos v ( a u+cu ) =vf (31) aplicando el teorema de Green (79) obtenemos la Ec. (??), que podemos reescribir como ( ) v a u+cuv dx= vf dx. (32) 1617 Definiendo el operador bilineal ( ) a(u,v)= v a u+cuv dx (33) y la funcional lineal l(v)= f,v = vf dx (34) podemosreescribirelproblemadadoporlaec.(??)deorden2,haciendouso delaformabilineala(, )ylafuncionallineall( ). Entonces entenderemos en el presente contexto un problema variacional con valores de frontera (VBVP) por uno de la forma: hallar una función u que pertenezcaaunespaciodehilbertv =H 1 0()yquesatisfagalaecuación a(u,v)= f,v (35) paratodafunciónv V dondea(, )esunaformabilinealyl( )esunafuncional lineal. Definición24 Sea V un espacio de Hilbert y sea V la norma asociada a dicho espacio, decimos que una forma bilineal a(, ) es continua si existe una constantem>0talque a(u,v) M u V v V u,v V (36) yesv-elípticosiexisteunaconstanteα>0talque donde V eslanormaasociadaalespaciov. a(v,v) α v 2 V v V (37) EstosignificaqueunaformaV elípticoesunaquesiempreesnonegativa ytomaelvalorde0sóloenelcasodequev=0,i.e.espositivadefinida. Notemos que el problema (30) definido en V = H 1 0() reescrito como el problema(35) genera una forma bilineal V-elíptico cuyo producto interior sobre V essimétricoypositivodefinidoyaque a(v,v) α v 2 V >0, v V,v 0 (38) reescribiéndoseelproblema(35),enelcualdebemosencontraru V talque dondeu 0 =gen,paratodav V. a(u,v)= f,v a(u 0,v) (39) Entonces, la cuestión fundamental, es conocer bajo que condiciones el problema anterior tiene solución y esta es única, el teorema de Lax-Milgram nos da las condiciones bajo las cuales el problema(30) reescrito como el problema(35) tiene solución y esta es única, esto queda plasmado en el siguiente resultado. 1718 Teorema 25 (Lax-Milgram) Sea V un espacio de Hilbert y sea a(, ) : V V R una forma bilineal continua V-elíptico sobre V. Además, sea l( ) : V R una funcional lineal continua sobre V. Entonces i)elvbvpdeencontraru V quesatisfaga a(u,v)= f,v, v V (40) tieneunaysólounasolución; ii)lasolucióndependecontinuamentedelosdatos,enelsentidodeque u V 1 α l V (41) donde V eslanormaenelespaciodualv dev yαes laconstantedela definición de V-elíptico. Másespecíficamente,considerandoahoraV unsubespaciocerradodeh m () las condiciones para la existencia, unicidad y la dependencia continua de los datos queda de manifiesto en el siguiente resultado. Teorema26 SeaV unsubespaciocerradodeh m (),seaa(, ):V V R unaformabilinealcontinuav-elípticosobrev yseal( ):V Runafuncional linealcontinuasobrev.seap unsubespaciocerradodev talque a(u+p,v+p)=a(u,v) u,v V yp,p P. (42) También denotando por Q el subespacio de V consistente de las funciones ortogonalesap enlanormal 2 ;talque { } Q= v V updx=0 p P, (43) yasumiendoquea(, )esq-elíptico:existeunaconstanteα>0talque a(q,q) α q 2 Q paraq Q, (44) lanormasobreqserálamismaquesobrev.entonces i)existeunaúnicasoluciónalproblemadeencontraru Qtalque si y sólo si las condiciones de compatibilidad se satisfacen. ii) La solución u satisface (dependencia continua de los datos). a(u,v)= l,v, v V (45) l,p =0 parap P (46) u Q α 1 l Q (47) 1819 Otroaspectoimportanteeslaregularidaddelasolución,silasoluciónual VBVPdeorden2mconf H s 2m ()dondes 2m,entoncesupertenecerá ah s ()yestoquedademanifiestoenelsiguienteresultado. Teorema27 Sea R n undominiosuaveyseau V lasoluciónalvbvp a(u,v)= f,v, v V (48) donde V H m (). Si f H s 2m () con s 2m, entonces u H s () y la estimación u H s C f H s 2m (49) se satisface. 1920 3. Solución de Grandes Sistemas de Ecuaciones Es este trabajo se muestra como proceder para transformar un problema de ecuaciones diferenciales parciales en un sistema virtual algebraico de ecuaciones lineales y así poder hallar la solución resolviendo el sistema lineal asociado al DDM, en general se pueden expresar en la forma matricial siguiente Au=f (50) dondelamatrizaesdetamañon nygeneralmentebandada,cuyotamaño debandaesb. En general, los métodos de resolución del sistema algebraico de ecuaciones Au=f seclasificanendosgrandesgrupos:losmétodosdirectosylosmétodos iterativos.enlosmétodosdirectoslasoluciónuseobtieneenunnúmerofijode pasos y sólo están sujetos a los errores de redondeo. En los métodos iterativos, se realizan iteraciones para aproximarse a la solución u aprovechando las característicaspropiasdelamatriza,tratandodeusarunmenornúmerodepasos queenunmétododirecto[7]a[11]. Para poder operar con los diversos métodos numéricos directos o iterativos que resuelven sistemas lineales, se programo una clase abstracta Solvable a la cual los métodos DDM le solicitan que resuelva un determinado sistema lineal sin importar el método numérico subyacente, la clase Solvable es una clase abstracta o contenedor con la siguiente estructura: class Solvable { private: int iter; public: virtual int getiter(void)=0; virtual void solve(double*x, double*y)=0; }; Lo cual permite operar con sistemas lineales virtuales o reales, usando la implementación de diversos métodos numéricos directos o iterativos, tanto secuenciales como paralelos indistintamente de las características propias de cada uno de ellos Métodos Directos Enlosmétodosdirectos,lasoluciónuseobtieneenunnúmerofijodepasos y sólo están sujetos a errores de redondeo. Entre los métodos más importantes podemos encontrar: Factorización LU-para matrices simétricas y no simétricasy descomposición de Cholesky-para matrices simétricas-. En todos los casos la matrizoriginalaesmodificadayencasodeusarlafactorizaciónlueltamaño delabandabcrecea2b+1silafactorizaciónserealizaenlamismamatriz. 2021 Descomposición LU Sea U una matriz triangular superior obtenida de A por eliminación bandada.entoncesu =L 1 A,dondeLesunamatriztriangularinferiorconunos enladiagonal.lasentradasdel 1 puedenobtenersedeloscoeficientesl ij y pueden ser almacenados estrictamente en las entradas de la diagonal inferior de A ya que estas ya fueron eliminadas. Esto proporciona una Factorización LU deaenlamismamatrizaahorrandoespaciodememoria,dondeelanchode bandacambiadeba2b+1. En el algoritmo de Factorización LU, se toma como datos de entrada del sistema Au=f (51) a la matriz A la cual será factorizada en una sola matriz la cual contendrá a las matrices L y U producto de la factorización, quedando el método numérico esquemáticamente como: A ii =1, parai=1,...,n ParaJ=1,2,...,N { Parai=1,2,...,j i 1 A ij =A ij A ik A kj } k=1 Parai=j+1,j+2,...,N ( ) j 1 A ij = 1 A jj A ij A ik A kj k=1 (52) ElproblemaoriginalAu=fseescribecomoLUu=f,dondelabusquedade la solución u se reduce a la solución sucesiva de los sistemas lineales triangulares i,.e. Ly=f y Uu=y. (53) y Ly = f y i = 1 A ii y 1 =f 1 /l 11 i 1 f i j=1 A ii y i paratodai=1,...,n (54) Uu = y x i = 1 A ii y i n j=i+1 x n =y n /u nn A ij x j paratodai=n 1,...,1.(55) 2122 LadescomposiciónLU requieren 3 /3operacionesaritméticasparalamatriz llena, pero sólo Nb 2 operaciones aritméticas para la matriz con un ancho de banda de b siendo esto más económico computacionalmente. Implementación LU Orientada a Objetos en C++ En este caso desarrolle una jerarquía de clases para poder operar con matrices bandadas y dispersas, llamada BandSolve la cual se muestra a continuación: class BandSolve: public Solvable { private: int bw; intn; double**ak; protect: factorlu(void); public: BanSolve(int n, double**a); BanSolve(int n, MatrizDispersa*A); void solve(double*x, double*y); void convertband(void); int getiter(void); }; Descomposición de Cholesky Cuando la matriz es simétrica y definida positiva, se obtiene la descomposición LU de la matriz A=LDU =LDL T donde D =diag(u) es la diagonal con entradas positivas. En el algoritmo de descomposición de Cholesky, se toma como datos de entrada del sistema Au=f (56) ala matriz A lacual seráfactorizada en unasolamatriz, lacual contendráa lasmatriceslyl T productodelafactorización,quedandoelmétodonumérico esquemáticamente como: parai=1,2,...,nyj=i+1,...,n ( ) i 1 A ii = A ii k=1 i 1 A 2 ik ( A ji = A ji A jk A ik )/A ii k=1 (57) 2223 ElproblemaoriginalAu=f seescribecomoll T u=b,dondelabusqueda delasoluciónusereducealasoluciónsucesivadelossistemaslinealestriangulares Ly=f y L T u=y (58) usando la formulación equivalente dada por las Ec.(54) y(55) para la descomposición LU. Lamayorventajadeestadescomposiciónesque,enelcasoenqueesaplicable,elcostodecómputoessustancialmentereducido,yaquerequieredeN 3 /6 multiplicacionesyn 3 /6adiciones. Implementación Cholesky Orientada a Objetos en C++ En este caso desarrolle una jerarquía de clases para poder operar con matrices bandadas y dispersas, llamada BandCholesky la cual se muestra a continuación: class BandCholesky: public Solvable { private: int bw; intn; double**ak; protect: factorlu(void); public: BanCholesky(int n, double**a); BanCholesky(int n, MatrizDispersa*A); void solve(double*x, double*y); void convertband(void); int getiter(void); }; 3.2. Métodos Iterativos En los métodos iterativos, se realizan iteraciones para aproximarse a la solución u aprovechando las características propias de la matriz A, tratando de usar unmenornúmerodepasosqueenunmétododirecto[7]y[10]. En los métodos iterativos tales como Jacobi, Gauss-Seidel y de Relajación Sucesiva(SOR) en el cual se resuelve el sistema lineal Au=f (59) comienzaconunaaproximacióninicialu 0 alasoluciónuygeneraunasucesión de vectores { u k} que converge a u. Los métodos iterativos traen consigo k=1 un proceso que convierte el sistema Au = f en otro equivalente mediante la iteracióndepuntofijodelaformau=tu+cparaalgunamatrizfijat yun 2324 vectorc.luegodeseleccionarelvectorinicialu 0 lasucesióndelosvectoresde la solución aproximada se genera calculando u k =Tu k 1 +c k=1,2,3,... (60) La convergencia a la solución la garantiza el siguiente teorema cuya solución puede verse en[11]. Teorema28 Si T <1,entonceselsistemalinealu=Tu+ctieneunasoluciónúnicau ylasiteracionesu k definidasporlafórmulau k =Tu k 1 +c k= 1,2,3,...convergenhacialasoluciónexactau paracualquieraproximaciónlinealu 0. Note que, mientras menor sea la norma de la matriz T, más rápida es la convergencia,enelcasocuando T esmenorqueuno,perocercanoauno,la convergencia es muy lenta y el número de iteraciones necesario para disminuir el error depende significativamente del error inicial. En este caso, es deseable proponer al vector inicial u 0 de forma tal que se mínimo el error inicial. Sin embargo,laeleccióndedichovectornotieneimportanciasila T espequeña, ya que la convergencia es rápida. Como es conocido, la velocidad de convergencia de los métodos iterativos dependen de las propiedades espectrales de la matriz de coeficientes del sistema de ecuaciones, cuando el operador diferencial L de la ecuación del problema a resolver es auto-adjunto se obtiene una matriz simétrica y positivo definida y el número de condicionamiento de la matriz A, es por definición cond(a)= λ máx λ mín 1 (61) donde λ máx y λ mín es el máximo y mínimo de los eigen-valores de la matriz A. Si el númerode condicionamientoes cercanoa1losmétodos numéricosal solucionar el problema convergerá en pocas iteraciones, en caso contrario se requerirán muchas iteraciones. Frecuentemente al usar el método de elemento finitosetieneunavelocidaddeconvergenciadeo ( ) 1 h yenelcasodemétodos 2 dedescomposicióndedominiosetieneunavelocidaddeconvergenciadeo ( ) 1 h en el mejor de los casos, donde h es la máxima distancia de separación entre nodos continuos de la partición, es decir, que poseen una pobre velocidad de convergenciacuandoh 0,paramásdetallesver[5],[6],[7] y[11]. Los métodos Jacobi, Gauss-Seidel y de Relajación Sucesiva(SOR) son usualmente menos eficientes que los métodos discutidos en el resto de esta sección basados en el espacio de Krylov[40]. Estos métodos minimizan, en la k-ésima iteraciónalgunamedidadeerrorsobreelespacioafinx 0 +K k,dondex 0 esla iteracióninicialyk k eselk-ésimosubespaciodekrylov } K k =Generado {r 0,Ar 0,...,A k 1 r 0 parak 1. Elresidualesr=b Ax,tal { r k }k 0 denotalasucesiónderesiduales r k =b Ax k. 2425 EntrelosmétodosmásusadosdefinidosenelespaciodeKrylovparaeltipo de problemas tratados en el presente trabajo podemos encontrar: Método de Gradiente Conjugado -para matrices simétricas- y GMRES -para matrices no simétricas Método de Gradiente Conjugado Silamatrizgeneradaporladiscretizaciónessimétrica ( A=A T) ydefinida positiva ( u T Au>0paratodou 0 ),entoncesesaplicableelmétododegradiente Conjugado (Conjugate Gradient Method CGM). La idea básica en que descansa el método del Gradiente Conjugado consiste en construir una base de vectores ortogonales espacio de Krylov K n ( A,v n ) y utilizarla para realizar la búsqueda de la solución en forma lo más eficiente posible. Tal forma de proceder generalmente no sería aconsejable porqué la construcción de una base ortogonal utilizando el procedimiento de Gramm-Schmidt requiere, al seleccionar cada nuevo elemento de la base, asegurar su ortogonalidad con respecto a cada uno de los vectores construidos previamente. La gran ventaja del método de Gradiente Conjugado radica en que cuando se utiliza este procedimiento, basta con asegurarlaortogonalidaddeunnuevomiembroconrespectoalúltimoqueseha construido, para que automáticamente esta condición se cumpla con respecto a todos los anteriores. EnelalgoritmodeGradienteConjugado,setomaalamatrizAcomosimétricaypositivadefinida,ycomodatosdeentradadelsistema Au=f (62) elvectordebúsquedainicialu 0 ysecalcular 0 =f Au 0,p 0 =r 0,quedando el método numérico esquemáticamente como: α n = pn,p n p n,ap n u n+1 =u n +α n p n r n+1 =r n α n Ap n Prueba de convergencia β n = rn+1,r n+1 r n,r n p n+1 =r n+1 +β n p n n=n+1 (63) donde, =(, )seráel productointerioradecuadoalsistemalinealenparticular,lasoluciónaproximadaseráu n+1 yelvectorresidualserár n+1. En la implementación numérica y computacional del método es necesario realizar la menor cantidad de operaciones posibles por iteración, en particular enap n,unamaneradehacerloquedaesquemáticamentecomo: Dadoelvectordebúsquedainicialu,calcular f Au,p ryµ r r. Mientras(µ<ε){ v Ap 2526 } α µ p v u u+αp r r αv µ r r β µ µ p r+βp µ µ n=n+1 Lasoluciónaproximadaseráuyelvectorresidualserár. Si denotamos {λ i,v i } N i=1 como las eigen-soluciones de A, i.e. AV i = λ i V i, i=0,1,2,...,n.yaquelamatrizaessimétrica,loseigen-valoressonrealesy podemosordenarlosporλ 1 λ 2... λ N.Definimoselnúmerodecondición porcond(a)=λ N /λ 1 ylanormadelaenergíaasociadaaapor u 2 A =u Au entonces u u k A u u 0 1 2k Cond(A) A. (64) 1+ Cond(A) El siguiente teorema nos da idea del espectro de convergencia del sistema Au=bparaelmétododeGradienteConjugado. Teorema29 Sea κ=cond(a)= λ máx λ mín 1, entonces el método de Gradiente Conjugado satisface la A norma del error dado por e n ( ) n e 0 2 κ 1 [ ( κ+1 ) n+ ( κ+1 ) ] n 2 (65) κ+1 κ 1 κ 1 dondee m =u u m delsistemaau=b. Notemosqueparaκgrandesetieneque κ 1 κ κ (66) tal que e n A e 0 A exp ( 2 κ n ) (67) deloanteriorpodemosesperarunespectrodeconvergenciadelordendeo( κ) iteraciones, para mayor referencia ver[11] y[40]. Definición30 Unmétodoiterativoparalasolucióndeunsistemalinealesllamado óptimo, si la razón de convergencia a la solución exacta es independiente del tamaño del sistema lineal. 2627 Definición 31 Un método para la solución del sistema lineal generado por métodos de descomposición de dominio es llamado escalable, si la razón de convergencia no se deteriora cuando el número de subdominios crece. La gran ventaja de los métodos de descomposición de dominio precondicionados entre los quedestacanafeti-dp ybddcesque estos métodos pueden ser óptimos y escalables según el tipo de precondicionador a priori que se use en ellos. ConsideracionessobrelaMatrizyelproductopuntoenelCGM Es común al trabajar en métodos de descomposición de dominio que se requieran usar más de un producto interior, en el mismo algoritmo de Gradiente Conjugado. Por ello se diseña el algoritmo para soportar en forma de parámetro el objeto que implemente el producto punto. class DotProd { public: virtual double dot(double*x, double*y)=0; }; EnnuestrocasoalusarDDMlamatrizconlaquetrabajaraelmétodode GradienteConjugadopuedeestardefinidacomotaloservirtual,oenelcaso más general, puede estar en un procesador o en múltiples procesadores, por ello en el diseño del algoritmo se deberá de implementar un objeto en el cual se pueda realizar la multiplicación de la matriz por el vector, sin importar el tipo de matriz-real o virtual-. class MultOp { public: virtual void multop(double*x, double*y)=0; virtual int getsize(void)=0; }; Así, se diseña un solo método de Gradiente Conjugado robusto y flexible que soporta diferentes productos interiores y trabaja con matrices reales o virtuales, adaptándose sin ningún cambio a diversas necesidades de implementación tanto secuencial como paralela. Implementación CGM Orientada a Objetos en C++ Dado el sistema linealau=f,dondeaseaunamatrizarealovirtual.laentradaalmétodoserá unaeleccióndeu 0 comocondicióninicial,ε>0comolatoleranciadelmétodo, N comoelnúmeromáximodeiteracionesademásdeunobjetopararealizarla multiplicación de la matriz por un vector y otro objeto para el producto interior, el algoritmo del método de Gradiente Conjugado orientado a objetos en C++ queda como: 2728 class CGM: public Solvable { protect: intn; double*r,*p,*v; double alpha, beta, gamma, lambda; double eps, mu; MultOp*A; DotProd*dotP; int nmaxiter; double norm(double*x); public: CGM(MulOp&A, DotProd&dotP, double eps); void solve(double*x, double*y); int getiter(void); }; Método Residual Mínimo Generalizado Si la matriz generada por la discretización es no simétrica, entonces es una opción es el método Residual Mínimo Generalizado (Generalized Minimum Residual Method GMRES), este representa una formulación iterativa común satisfaciendo una condición de optimización. La idea básica detrás del método se basa en construir una base ortonormal { v 1,v 2,...,v n} (68) ( paraelespaciodekrylovk ) ( n A,v n.parahacerv n+1 ortogonalak ) n A,v n, esnecesariousartodoslosvectorespreviamenteconstruidos { v n+1j} n j=1 (enla práctica sólo se guardan unos cuantos vectores anteriores) en los cálculos. Y el algoritmo se basa en una modificación del método de Gram-Schmidt para la generacióndeunabaseortonormal.seav n = [ v 1,v 2,...,v n] lacualdenotala matriz conteniendov j enlaj-ésimacolumna,paraj=1,2,...,n,yseah n = [h i,j ],1 i,j n,donde lasentradasde H n noespecificadasenel algoritmo soncero.entonces,h n esunamatrizsuperiordehessenberg.i.e.h ij =0para j<i 1,y AV n = V n H n +h n+1,n [ 0,...,0,v n+1 ] (69) H n = H T n AV n En el algoritmo del método Residual Mínimo Generalizado, la matriz A es tomada como no simétrica, y como datos de entrada del sistema Au=f (70) 2829 el vector de búsqueda inicial u 0 y se calcula r 0 = f Au 0, β 0 = r 0, v 1 = r 0 /β 0,quedandoelmétodoesquemáticamentecomo: Paran=1,2,...,Mientrasβ n <τβ 0 { w n+1 0 =Av n Paral=1hastan{ } h l,n = w n+1 l,v l w n+1 l+1 =wn+1 l h l,n v l } h n+1,n = w n+1 n+1 v n+1 =w n+1 n+1 /h n+1,n Calculary n talqueβ n = β 0 e 1 Ĥn yn esmínima (71) dondeĥn =[h ij] 1 i n+1,1 j n,lasoluciónaproximadaseráu n =u 0 +V n y n, y el vector residual será r n =r 0 AV n y n =V n+1 (β 0 e 1 Ĥ n y n). Teorema32 Seau k laiteracióngeneradadespuésdekiteracionesdegmres, con residual r k. Silamatriz Aes diagonalizable, i.e.a=vλv 1 dondeλes unamatrizdiagonaldeeigen-valoresdea,yv eslamatrizcuyascolumnasson los eigen-vectores, entonces r k r 0 κ(v) mín máx p κ(λ j ) p κ Πκ,p k (0)=1 λ j dondeκ(v)= V V 1 eselnúmerodecondicionamientodev. Consideraciones sobre la Matriz en el GMRES En nuestro caso al usar DDM la matriz con la que trabajara el método GMRES puede estar definida comotaloservirtual,oenelcasomásgeneral,puedeestarenunprocesador o en múltiples procesadores, por ello en el diseño del algoritmo se deberá de implementar un objeto en el cual se pueda realizar la multiplicación de la matriz porelvector,sinimportareltipodematriz-realovirtual-. class MultOp { public: virtual void multop(double*x, double*y)=0; virtual int getsize(void)=0; }; Así,sediseñaunsólométododeGMRESrobustoyflexiblequetrabajecon matrices reales o virtuales, adaptándose sin ningún cambio a diversas necesidades de implementación tanto secuencial como paralela. 2930 Implementación GMRES Orientada a Objetos en C++ Dado el sistema lineal Au = f, donde A sea una matriza real o virtual. La entrada al métodoseráunaeleccióndeu 0 comocondicióninicial,ε>0comolatolerancia delmétodo,ncomoelnúmeromáximodeiteracionesykelnúmerodevectores aguardarademásdeunobjetopararealizarlamultiplicacióndelamatrizpor un vector, el algoritmo del método de GMRES orientado a objetos en C++ queda como: class DQGMRES: public Solvable { protect: intn; intk; int MaxIter; MultOp*A; double gm, gm1; double**p; double**cs; double**h; double**q; double*v; double eps; public: DQGMRES(MulOp&A, int k); void void solve(double*x, double*y); int getiter(void); }; 3.3. Estructura óptima de las Matrices en su Implementación Computacional Una parte importante de la implementación computacional de los métodos numéricos de resolución de sistemas algebraicos en C++ es utilizar una forma óptima de almacenar, recuperar y operar las matrices, tal que, facilite los cálculos que involucra la resolución de grandes sistemas de ecuaciones cuya implementación puede ser secuencial o paralela. El sistema lineal puede ser expresado en la forma matricial Au=f (72) dondelamatriza(quepuedeserrealovirtual)esdetamañon nconbandab, peroelnúmerototaldedatosalmacenadosenellaesalosmásn bnúmerosde dobleprecisión,enelcasodesersimétricalamatriz,elnúmerodedatosalmacenadosesmenora(n b)/2.ademássielproblemaquelaoriginóesdecoeficientes constantes el número de valores almacenados se reduce drásticamente a sólo el tamañodelabandab. 3031 En el caso de que el método para la resolución del sistema lineal a usar sea del tipo Factorización LU o Cholesky, la estructura de la matriz cambia, ampliándoseeltamañodelabandadeba2 b+1enlafactorización,enelcaso de usar métodos iterativos tipo CGM o GMRES la matriz se mantiene intacta conunabandab. Para la resolución de los métodos de descomposición de dominio en los cuales es necesario usar métodos numéricos algebraicos para resolver los sistemas lineales asociados (y en este proceso se consume casi la totalidad del tiempo de ejecución), la operación básica que se realiza de manera reiterada es la multiplicacióndeunamatrizporunvectorv=cu,delaformamáseficienteposible. Un factor determinante en la eficiente implementación computacional es la forma de almacenar, recuperar y realizar las operaciones que involucren matrices yvectores,detalformaquelamultiplicaciónserealiceenlamenorcantidadde operaciones y que los valores necesarios para realizar dichas operaciones queden en la medida de lo posible contiguos para ser almacenados en el Cache del procesador. Dado que la multiplicación de una matriz C por un vector u, dejando el resultado en v se realiza mediante el algoritmo for(i=0; i<ren; i++) { s=0.0; for(j=0;j<col;j++) { s+=c[i][j]*u[j]; } v[i]=s; } Es necesario para lograr una eficiente implementación del algoritmo anterior, queelgranvolumendedatosdesplazadosdelamemoriaalcacheyviceversa durante la ejecución del programa sea óptimo. Notemos que la velocidad de acceso a la memoria principal es relativamente lento con respecto al Cache, el cual está dividido en L1, L2 y hasta L3, siendo el más rápido el L1 y el más lento el L3. Así, cada vez que las unidades funcionales requieren un conjunto de datos para implementar una determinada operación en los registros lo solicitan primeramente a los Caches, estos consumen diversa cantidad de ciclos de reloj paraentregareldatosilotienen,perosiempreeltiempoesmenorquesolicitarle el dato a la memoria principal. Por ello, los datos se deben agrupar para que la operación más usada -en nuestro caso multiplicación matriz vector- se realice con la menor solicitud de datosalamemoriaprincipal,silosdatosusados-renglóndelamatriz-seponen contiguos minimizará los accesos a la memoria principal, pues es más probable que estos estarán contiguos en el Cache al momento de realizar la multiplicación. Por ejemplo, en el caso de matrices bandadas de tamaño de banda b, el algoritmo anterior se simplifica a 3132 for(i=0; i<ren; i++) { s= 0.0; for(k=0;k<ban;k++) { if((ind[k]+i)>=0&&(ind[k]+i)<ren) s+=dat[i][k]*u[ind[k]+i]; } v[i]=s; } Si,lasolicituddememoriaparaDat[i]sehacedetalformaquelosdatosdel renglón estén continuos(son b números de punto flotante), esto minimizará los accesos a la memoria principal en cada una de las operaciones involucradas en el producto, como se explica en las siguientes secciones Matrices Bandadas Las matrices bandadas de banda b-sin perdida de generalidad y para propósitos de ejemplificación la supondremos pentadiagonal- surgida de descomposición,esdeltipo A= a 1 b 1 c 1 d 2 a 2 b 2 c 2 d 3 a 3 b 3 c 3 d 4 a 4 b 4 c 4 e 5 d 5 a 5 b 5 c 5 e 6 d 6 a 6 b 6 e 7 d 7 a 7 b 7 e 8 d 8 a 8 b 8 e 9 d 9 a 9 la cual puede ser almacenada usando el algoritmo Compressed Diagonal Storage (CDS)[10], optimizado para ser usado en C a 1 b 1 c 1 0 d 2 a 2 b 2 c 2 0 d 3 a 3 b 3 c 3 0 d 4 a 4 b 4 c 4 Dat= e 5 d 5 a 5 b 5 c 5 e 6 d 6 a 6 b 6 0 e 7 d 7 a 7 b 7 0 e 8 d 8 a 8 b 8 0 e 9 d 9 a conunvectordeíndices Ind=[ 5, 1,0,+1,+5] 3233 detalformaquepuedeserreconstruidalamatrizadeformaeficiente. ParaobtenerelvalorA i,j,calculamosind=j i,sielvalorindestaenla lista de índices Ind, -supóngase en la columna k entonces A i,j = Dat ik, en otrocasoa i,j =0. Casos Particulares de la Matriz Bandada A Básicamente dos casos particulares surgen en el tratamiento de ecuaciones diferenciales parciales: El primer casoes cuando el operadordiferencial parcial es simétricoyel otro, enel que los coeficientes del operador sean constantes. Para el primer caso, al ser la matriz simétrica, sólo es necesario almacenar la parte con índices mayores o iguales a cero, de tal forma que se buscara el índicequesatisfagaind= j i,reduciendoeltamañodelabandaab/2enla matriz A. Para el segundo caso, al tener coeficientes constantes el operador diferencial, losvaloresdelosrenglonesdentrodecadacolumnadelamatrizsoniguales,y sólo es necesario almacenarlos una sola vez, reduciendo drásticamente el tamaño delamatrizdedatos. Implementación de Matrices Bandadas Orientada a Objetos en C++ UnaformadeimplementarlamatrizbandadaAdebandabenC++esmediante: //Creacion double**dat; Dat=newdouble*[ren]; for(i=0;i<ren;i++)dat[i]=newdouble[b]; int*ind; Ind=newint[b]; // Inicializacion for(i=0;i<ren;i++) for(j=0;j<b;j++)dat[i][j]=0.0; for(i=0;i<b;i++)ind[i]=0; De esta forma, al tomar en cuenta la operación de multiplicación de una matriz por un vector, donde el renglón de la matriz involucrado en la multiplicación queda generalmente en una región contigua del Cache, se hace óptima la multiplicación de la matriz por el vector Matrices Dispersas Las matrices dispersas de a lo más b valores distintos por renglón-supongamos b=3-quesurgeenmétodosdedescomposicióndedominioparaalmacenaral- 3334 gunasmatricesysondeltipo a 1 b 1 c 1 a 2 b 2 c 2 a 3 b 3 c 3 a 4 b 4 A= a 5 b 5 c 5 a 6 b 6 c 6 a 7 b 7 c 7 a 8 b 8 c 8 a 9 b 9 la cual puede ser almacenada usando Jagged Diagonal Storage(JDC)[10], optimizadoparaserusadoenc++ a 1 b 1 c 1 a 2 b 2 c 2 a 3 b 3 c 3 a 4 b 4 0 Dat= a 5 b 5 c 5 a 6 b 6 c 6 a 7 b 7 c 7 a 8 b 8 c 8 a 9 b 9 0 con una matriz de índices Ind= detalformaquepuedeserreconstruidalamatrizadeformaeficiente. Para obtener el obtener el valor A i,j, buscamos el valor j en la lista de índices Ind dentro del renglón i, si lo encontramos en la posición k, entonces A i,j =Dat ik,enotrocasoa i,j =0. CasosParticularesdelaMatrizDispersaA SilamatrizA,quealserse almacenada, se observa que existen a lo más r diferentes renglones con valores distintosdelosnconquecuentalamatrizysir<<n,entoncesesposiblesólo guardar los r renglones distintos y llevar un arreglo que contenga la referencia al renglón almacenado. 3435 Implementación de Matrices Dispersas Orientada a Objetos en C++ UnaformadeimplementarlamatrizbandadaAdebandabenC++esmediante: //Creación double**dat; Dat=newdouble*[ren]; for(i=0;i<ren;i++)dat[i]=newdouble[b]; int**ind; Ind=newint*[ren]; for(i=0;i<ren;i++)ind[i]=newint[b]; // Inicialización for(i=0;i<ren;i++) for(j=0;j<b;j++)dat[i][j]=0.0,ind[i][j]=-1; De esta forma, al tomar en cuenta la operación de multiplicación de una matriz por un vector, donde el renglón de la matriz involucrado en la multiplicación queda generalmente en una región contigua del Cache, se hace óptima la multiplicación de la matriz por el vector Multiplicación Matriz-Vector Los métodos de descomposición de dominio requieren por un lado la resolución de almenos un sistema lineal y por el otro lado requieren realizar la operacióndemultiplicacióndematriz porvector,i.e. Cudelaformamáseficiente posible, por ello los datos se almacenan de tal forma que la multiplicación se realice en la menor cantidad de operaciones. Dado que la multiplicación de una matriz C por un vector u, dejando el resultado en v se realiza mediante el algoritmo: for(i=0; i<ren; i++) { s=0.0; for(j=0;j<col;j++) { s+=c[i][j]*u[j]; } v[i]=s; } Enelcasodematricesbandadas,sesimplificaa: for(i=0; i<ren; i++) { s= 0.0; for(k=0;k<ban;k++) { if((ind[k]+i)>=0&&(ind[k]+i)<ren) 3536 } } v[i]=s; s+=dat[i][k]*u[ind[k]+i]; De forma similar, en el caso de matrices dispersas, se simplifica a: for(i=0; i<ren; i++) { s=0.0,k=0 while(ind[i][k]!=-1) { s+=dat[i][k]*u[ind[i][k]]; k++; if(k>=b)break; } v[i]=s; } De esta forma, al tomar en cuenta la operación de multiplicación de una matriz por un vector, donde el renglón de la matriz involucrado en la multiplicación queda generalmente en una región contigua del Cache se hace óptima la multiplicación de la matriz por el vector. 3637 4. Métodos de Solución Aproximada para EDP En el presente capítulo se prestará atención a varios aspectos necesarios para encontrar la solución aproximada de problemas variacionales con valor en la frontera (VBVP). Ya que en general encontrar la solución a problemas con geometría diversa es difícil y en algunos casos imposible usando métodos analíticos. EnestecapítuloseconsideraelVBVPdelaforma donde Lu = f en (73) u = g en Lu= a u+cu (74) conaunamatrizpositivadefinida,simétricayc 0,comouncasoparticular deloperadorelípticodefinidoporlaec.(??)deorden2,con R 2 undominio poligonal,esdecir,esunconjuntoabiertoacotadoyconexotalquesufrontera eslaunióndeunnúmerofinitodepolígonos. La sencillez del operador L nos permite facilitar la comprensión de muchas de las ideas básicas que se expondrán a continuación, pero tengamos en mente que esta es una ecuación que gobierna los modelos de muchos sistemas de la ciencia y la ingeniería, por ello es muy importante su solución. Simultiplicamosalaecuación a u+cu=f porv V =H 1 0(), obtenemos v ( a u+cu ) =vf (75) aplicando el teorema de Green (79) obtenemos la Ec. (??), que podemos reescribir como ( ) v a u+cuv dx= vf dx. (76) Definiendo el operador bilineal ( ) a(u,v)= v a u+cuv dx (77) y la funcional lineal l(v)= f,v = vf dx (78) podemosreescribirelproblemadadoporlaec.(73)deorden2enformavariacional,haciendousodelaformabilineala(, )ylafuncionallineall( ) Método Galerkin La idea básica detrás del método Galerkin es, considerando el VBVP, encontraru V =H 1 0 ()quesatisfaga a(u,v)= f,v v V (79) 3738 donde V es un subespacio de un espacio de Hilbert H (por conveniencia nos restringiremos a espacios definidos sobre los números reales). ElproblemaaltratarderesolverlaEc.(79)estáenelhechodequeelespacio V es de dimensión infinita, por lo que resulta que en general no es posible encontrarelconjuntosolución.enlugardetenerelproblemaenelespaciov, se supone que se tienen funciones linealmente independientes φ 1,φ 2,...,φ N en V ydefinimoselespaciov h apartirdelsubespaciodimensionalmentefinitode V generadoporlasfuncionesφ i,esdecir, V h =Generado{φ i } N i=1, Vh V. (80) Elíndiceh=1/N esunparámetroqueestaráentre0y1,cuyamagnitudda algunaindicacióndecuancercav h estadev,hserelacionaconladimensión dev h.ycomoelnúmeron delasfuncionesbaseseescogedemaneraquesea grandeyhagaquehseapequeño,enellímite,cuandon,h 0. DespuésdedefinirelespacioV h,esposibletrabajarconv h enlugardev yencontrarunafunciónu h quesatisfaga a(u h,v h )= f,v h v h V h. (81) EstaeslaesenciadelmétodoGalerkin,notemosqueu h yv h sonsólocombinacioneslinealesdelasfuncionesbasedev h,talesque N N u h = c i φ i y v h = d j φ j (82) i=1 dondev h esarbitraria,comoloscoeficientesded j ysinperdidadegeneralidad podemoshacerv h =φ j.así,paraencontrarlasoluciónu h sustituimoslasecs. (82)enlaEc.(81)yusandoelhechoquea(, )esunaformabilinealyl( )es una funcional lineal se obtiene la ecuación j=1 N a ( ) φ i,φ j ci = f,φ j i=1 (83) o más concisamente, como enlacual N K ij c i F j =0 j=1,2,...n (84) i=1 K ij =a ( φ i,φ j ) y F j = f,φ j (85) notemosquetantok ij yf j puedenserevaluados,yaqueφ i,a(, )yl( )son conocidas. Entonces el problema se reduce a resolver el sistema de ecuaciones lineales N K ij c i F j, j=1,2,...n (86) i=1 3839 o más compactamente Ku=F (87) enlacualkyf sonlamatrizyelvectorcuyasentradassonk ij yf j.unavez queelsistemaesresuelto,lasoluciónaproximadau h esencontrada. Notemos que la forma bilineal a(, ) define un producto interior sobre V, si a(, ) es simetrica y V elíptica, entonces las propiedades de linealidad y simetriasonobvias,mientrasquelapropiedaddev elíticidaddea(, )espor a(v,v) α v 2 >0 v 0, (88) además, si a(, ) es continua, entonces la norma v a a(v,v) generada por esteproductointerioresequivalentealanormaestandarsobrev,talquesiv es completa con respecto a la norma estandar, esta también es completa con respectoalanorma v a. Porotrolado,sielconjuntodefuncionesbase{φ i } N i=1 seeligendetalforma que sean sean ortogonales entre si, entonces el sistema (84) se simplifica considerablemente, ya que K ij =a ( φ i,φ j ) =0 sii j (89) y K ii c i =F i ó c i =F i /K ii. (90) Así,elproblema(73)definidoenV h =H 1 0 ()reescritocomoelproblema (79)generaunaformabilineal V h -elíptica cuyoproductointeriorsobre V h es simétrico y positivo definido ya que a(v h,v h ) α v h 2 V h >0, v h V h,v h 0 (91) reescribiéndose el problema(81) como el problema aproximado en el cual debemosencontraru h V h V talque a(u h,v h )= f,v h a(u 0,v h ) (92) dondeu 0 =g=0en,paratodav h V h,esdecir ( ) vh a u h +cu h v h dxdy= f v h dxdy (93) paratodov h V h. Entonces, el problema(73) al aplicarle el método Galerkin obtenemos(76), el cual podemos reescribirlo como(93). Aplicando el teorema de Lax-Milgram (25) a este caso particular, tenemos que este tiene solución única y esta depende continuamente de los datos. Como un caso particular del teorema de Lax-Milgram(25) tenemoe el siguiente resultado 3940 Teorema33 SeaV h unsubespaciodedimensiónfinitadeunespaciodehilbert V,seaa(, ):V h V h RunaformabilinealcontinuayV-elíptica,yl( ): V h Runafuncionallineal acotada.entoncesexisteunaúnicafuncionu h V h talquesatisface Además,sil( )esdelaforma l,v h = a(u h,v h )= l,v h v h V h. (94) f v h dx (95) conf L 2 (),entonces u h V 1 α f L2, (96) dondeαeslaconstanteen(88). El siguiente resultado nos da una condición suficiente para que la aproximaciónu h delmétodogalerkinconverjaalasoluciónudelproblemadadopor laec.(79),paramásdetallevéase[?]y[?]. Teorema34 Sea V un subespacio cerrado de un espacio de Hilbert, y sea la formabilineala(, ):V h V h RcontinuaV-elípticayseal( )unafuncional lineal acotada. Entonces existe una constante C, independiente de h, tal que u u h V C ínf v h V h u v h V (97) donde u es solución de (79) y u h es solución de (92), consecuentemente, una condiciónsuficienteparaquelaaproximaciónu h delmétodogalerkinconverge alasoluciónudelproblemadadoporlaec.(79)esqueexistaunafamilia { V h} de subespacios con la propiedad de que ínf u v h V 0 cuando h 0. (98) v h V h 4.2. El Método de Residuos Pesados EstemétodosebasaenelmétodoGalerkin,yseescogensubespaciosU h y V h detalmaneraqueladimensióndimu h =dimv h =N,eligiendolasbases como {φ i } N i=1 parauh y { ψ j } N j=1 paravh (99) entonces N N u h = c i φ i yv h = b j ψ j (100) i=1 dondeloscoeficientesb j sonarbitrariosyaquev h esarbitraria. j=1 4041 Sustituyendo esta última expresión Ec.(100) en se obtienen N ecuaciones simultaneas (Lu h f,v h )=0 (101) N K ij c i =F j conj=1,...,n i=1 enlacualenlacualkyf sonlamatrizyelvectorcuyasentradasson K ij = ( Lφ i,ψ j ) yfj = ( f,ψ j ) donde(, )representaelproductointeriorasociadoal 2.Alaexpresión τ(u h ) Lu h f (102) selellamaelresiduo;siu h eslasoluciónexacta,entoncesporsupuestoelresiduo se nulifica Método de Elementos Finitos El método Finite Elements Method(FEM) provee una manera sistemática y simple de generar las funciones base en un dominio con geometría poligonal. Lo que hace al método de elemento finito especialmente atractivo sobre otros métodos, es el hecho de que las funciones base son polinomios definidos por pedazos (elementos i ) que son no cero sólo en una pequeña parte de, proporcionando a la vez una gran ventaja computacional al método ya que las matrices generadas resultan bandadas ahorrando memoria al implantarlas en una computadora. Así, partiendo del problema aproximado(93), se elegirá una familia de espaciosv h (h (0,1))definidoporelprocedimientodeelementosfinitos(descritos en las subsecciones siguientes en el caso de interpoladores lineales, para otros tiposdeinterpoladores,ver[7]),teniendolapropiedaddequev h seaproxima av cuandohseaproximaaceroenunsentidoapropiado,estoes,porsupuesto una propiedad indispensable para la convergencia del método Galerkin. Mallado del dominio El Mallado o triangulación T h del dominio es el primeraspectobásico,yciertamenteelmáscaracterístico,eldominio R 2 essubdivididoenesubdominiosoelementos e llamadoselementosfinitos,tal que E = donde: e=1 e 4142 Cada e T h esunpolígono(rectángulootriángulo)coninterior novacío( e )yconexo. Cada e T h tienefrontera e Lipschitzcontinua. Paracada i, j T h distintos, i j =. Eldiámetroh i =Diam( e )decada e satisfacediam( e ) h paracadae=1,2,...,e. Cualquiercaradecualquierelemento i T h enlatriangulación es tambiénunsubconjuntode lafrontera del dominioouna caradecualquierotroelemento j T h delatriangulación,eneste últimocaso i y j sonllamadosadyacentes. Losvérticesdecada e sonllamadosnodos,teniendon deellos porcadaelemento e. Una vez que la triangulación T h del dominio es establecida, se procede a definir el espacio de elementos finitos P h [k] a través del proceso descrito a continuación. Funciones Base A continuación describiremos la manera de construir las funciones base usada por el método de elemento finito. En este procedimiento debemos tener en cuenta que las funciones base están definidas en un subespacio dev =H 1 ()paraproblemasdesegundoordenquesatisfacenlascondiciones de frontera. Las funciones base deberán satisfacer las siguientes propiedades: Lasfuncionesbaseφ i sonacotadasycontinuas,i.eφ i C( e ). Existen l funciones base por cada nodo del polígono e, y cada funciónφ i esnocerosoloenloselementoscontiguosconectadospor elnodoi. φ i =1encadainododelpolígono e yceroenlosotrosnodos. Larestricciónφ i a e esunpolinomio,i.e.φ i P k [ e ]paraalguna k 1dondeP k [ e ]eselespaciodepolinomiosdegradoalomásk sobre e. Decimos que φ i P k [ e ] es una base de funciones y por su construcción esevidentequeestaspertenecenah 1 ().Alconjuntoformadoportodaslas funciones base definidas para todo e de será el espacio P h [k] de funciones base, i.e. E P h [k]= P k [ e ] e=1 estas formarán las funciones base globales. 4243 Solución aproximada Para encontrar la solución aproximada elegimos el espaciop h [k]defuncionesbase,comoelespaciodefuncioneslinealesφ i definidas por pedazos de grado menor o igual a k (en nuestro caso k = 1), entonces el espacio a trabajar es V h =Generado { φ i P h [k] φ i (x)=0en }. (103) La solución aproximada de la Ec. (93) al problema dado por la Ec. (73) queda en términos de ( ) φi a φ j cφ i φ j dxdy= f φ j dxdy (104) si definimos el operador bilineal K ij a ( φ i,φ j ) = y la funcional lineal F j f,φ j = ( φi a ij φ j cφ i φ j ) dxdy (105) f φ j dxdy (106) entonceslamatrizk [K ij ],losvectoresu (u 1,...,u N )yf (F 1,...,F N ) definen el sistema lineal(que es positivo definido) Ku=F (107) dondeuseráelvectorsoluciónalaec.(107)cuyosvaloresseránlasoluciónal problemadadoporlaec.(93)que eslasoluciónaproximadaalaec. (73)en los nodos interiores de. Un Caso más General Sea el operador elíptico(caso simétrico) en el dominio,yeloperadordefinidopor donde Lu = f en\σ (108) u = g en [u] Σ = J 0 [a n u] Σ = J 1 Lu= a u+cu (109) conjuntamenteconunapartición ={ 1,..., E }de.multiplicandoporla función w obtenemos wlu= w a u+cwu=wf (110) entoncessiw(x)estalque[w]=0(esdecirwescontinua)ydefinimos 4344 a(u,w)= E i=1 talquea(u,w)defineunproductointeriorsobre i ( u a w+cwu ) dx (111) H 1 ()=H 1 ( 1 ) H 1 ( 2 )... H 1 ( E ). Entonces, reescribimos la Ec.(110) como a(u,w) = = wfdx+ E i=1 wf dx+ wa n uds (112) wa n uds Σ w[a n u] Σ ds. Sea u 0 (x) una función que satisface las condiciones de frontera y J 0 una función que satisface las condiciones de salto, tal que i)u 0 (x)=g(x)en ii)[u 0 (x)] Σ =J 0 yseau(x)=u 0 (x)+v(x).entoncesu(x)satisfacelaec.(111)siysólosi v(x) satisface a(u,w)= wf dx u 0,w J 1 wds (113) Σ paratoda w tal que w(x)=0en. Sea{φ i }una base de unsubespaciode dimensiónfinitav h definidocomo V h = { φ i φ i C 1 ( i ), i,φ i =0en yφ i C 0 () }. (114) La solución por elementos finitos de (113) se obtiene al resolver el sistema lineal Ku=F (115) donde y F j = K ij =a ( φ i,φ j ) φ j f dx a ( u 0,φ j ) Σ (116) J 1 φ j ds (117) esta solución será la solución en los nodos interiores de.en el presente capítulo se prestará atención a varios aspectos necesarios para encontrar la solución aproximada de problemas variacionales con valor en la frontera(vbvp). Ya que en general encontrar la solución a problemas con geometría diversa es difícil y en algunos casos imposible usando métodos analíticos. 4445 5. Método de Elementos Finitos donde EnestecapítuloseconsideraelVBVPdelaforma Lu = f en (118) u = g en Lu= a u+cu (119) conaunamatrizpositivadefinida,simétricayc 0,comouncasoparticular deloperadorelípticodefinidoporlaec.(??)deorden2,con R 2 undominio poligonal,esdecir,esunconjuntoabiertoacotadoyconexotalquesufrontera eslaunióndeunnúmerofinitodepolígonos. La sencillez del operador L nos permite facilitar la comprensión de muchas de las ideas básicas que se expondrán a continuación, pero tengamos en mente que esta es una ecuación que gobierna los modelos de muchos sistemas de la ciencia y la ingeniería, por ello es muy importante su solución. Simultiplicamosalaecuación a u+cu=f porv V =H 1 0 (), obtenemos v ( a u+cu ) =vf (120) aplicando el teorema de Green (79) obtenemos la Ec. (??), que podemos reescribir como ( ) v a u+cuv dx= vf dx. (121) Definiendo el operador bilineal ( ) a(u,v)= v a u+cuv dx (122) y la funcional lineal l(v)= f,v = vf dx (123) podemos reescribir el problema dado por la Ec. (118) de orden 2 en forma variacional,haciendousodelaformabilineala(, )ylafuncionallineall( ) Triangulación El Mallado o triangulación T h del dominio es el primer aspecto básico, y ciertamente el más característico, el dominio R 2 es subdividido en E subdominiosoelementos e llamadoselementosfinitos,talque donde: E = e=1 e 4546 Cada e T h esunpolígono(rectángulootriángulo)coninterior novacío( e )yconexo. Cada e T h tienefrontera e Lipschitzcontinua. Paracada i, j T h distintos, i j =. Eldiámetroh i =Diam( e )decada e satisfacediam( e ) h paracadae=1,2,...,e. Losvérticesdecada e sonllamadosnodos,teniendon deellos porcadaelemento e. Definición35 UnafamiliadetriangulacionesT h esllamadadeforma-regular si existe una constante independiente de h, tal que h K Cρ K, conk T h, dondeρ K es el radio del circulo más grandecontenido en K.El radio h K /ρ K esllamadoeslaspectratiodek. Definición36 UnafamiliadetriangulacionesT h esllamadacuasi-uniformesi estaesdeforma-regularysiexisteunaconstanteindependientedeh,talque h K Ch, conk T h. Una vez que la triangulación T h del dominio es establecida, se procede a definir el espacio de elementos finitos P h [k] a través del proceso descrito a continuación Interpolación para el Método de Elementos Finitos Funciones Base A continuación describiremos la manera de construir las funciones base usada por el método de elemento finito. En este procedimiento debemos tener en cuenta que las funciones base están definidas en un subespacio dev =H 1 ()paraproblemasdesegundoordenquesatisfacenlascondiciones de frontera. Las funciones base deberán satisfacer las siguientes propiedades: i)lasfuncionesbaseφ i sonacotadasycontinuas,i.eφ i C( e ). ii)existenlfuncionesbaseporcadanododelpolígono e,ycada funciónφ i esnocerosoloenloselementoscontiguosconectadospor elnodoi. iii)φ i =1encadainododelpolígono e yceroenlosotrosnodos. iv) La restricción φ i a e es un polinomio, i.e. φ i P k [ e ] para algunak 1dondeP k [ e ]eselespaciodepolinomiosdegradoalo másksobre e. 4647 Decimos que φ i P k [ e ] es una base de funciones y por su construcción esevidentequeestaspertenecenah 1 ().Alconjuntoformadoportodaslas funciones base definidas para todo e de será el espacio P h [k] de funciones base, i.e. E P h [k]= P k [ e ] e=1 estas formarán las funciones base globales Método de Elemento Finito Usando Discretización de Rectángulos Para resolver la Ec.(118), usando una discretización con rectángulos, primero dividimoseldominio R 2 enn x nodoshorizontalesporn y nodosverticales, teniendoe=(n x 1)(N y 1)subdominiosoelementosrectangulares e tales que= E e=1 e y i j sisonadyacentes,conuntotalden =N x N y nodos. Donde las funciones lineales definidas por pedazos en e en nuestro caso serán polinomios de orden uno en cada variable separadamente y cuya restricción deφ i a e esφ (e) i (. Para simplificar ) los cálculos en esta etapa, supondremos que 1 0 lamatriza=a,entoncessetienequelaintegraldelladoizquierdo 0 1 delaec.(174)quedaescritacomo ( ) a φi φ j +cφ i φ j dxdy= f φ j dxdy (124) donde K ij = = = E e=1 E e=1 ( a φi φ j +cφ i φ j ) dxdy (125) e e ( ) a φ (e) i φ (e) j +cφ (e) i φ (e) j dxdy ( a [ φ (e) φ (e) i j x + φ(e) i x y φ (e) ] ) j +cφ (e) i φ (e) j dxdy y yelladoderechocomo F j = = E e=1 f φ j dxdy (126) e f φ (e) j dxdy. 4748 Paracada e de,lasubmatrizdeintegrales(matrizdecargalocal) K ij = e ( tiene la estructura a [ φ (e) φ (e) i j x + φ(e) i x y K (e) 1,1 K (e) 1,2 K (e) 1,3 K (e) 1,4 K (e) 2,1 K (e) 2,2 K (e) 2,3 K (e) 2,4 K (e) 3,1 K (e) 3,2 K (e) 3,3 K (e) 3,4 K (e) 4,1 K (e) 4,2 K (e) 4,3 K (e) 4,4 φ (e) ] ) j +cφ (e) i φ (e) j dxdy (127) y lacualdeberáserensambladaenlamatrizdecargaglobalquecorrespondaa lanumeracióndenodos localesdel elemento e conrespecto alanumeración globaldeloselementosen. De manera parecida, para cada e de se genera el vector de integrales (vector de carga local) con la estructura F j = e f φ (e) j dxdy (128) F (e) 1 F (e) 2 F (e) 3 F (e) 4 el cual también deberá ser ensamblado en el vector de carga global que corresponda a la numeración de nodos locales al elemento e con respecto a la numeración global de los elementos de. Montando los K (e) ij en la matriz K y los F (e) j en el vector F según la numeracióndenodosglobal,segeneraelsistemaku h =Fdondeu h seráelvector cuyos valores serán la solución aproximada a la Ec.(118) en los nodos interiores de.lamatrizkgeneradadeestaforma,tieneunapropiedadmuyimportante, esbandadayelanchodebandaesde9elementos,estoesmuyútilalmomento de soportar la matriz en memoria. y Paraimplementarnuméricamenteencada e lasintegrales ( [ φ (e) φ (e) i j + φ(e) φ (e) ] ) i j a +cφ (e) i φ (e) j dxdy (129) x x y y e e f φ (e) j dxdy, (130) teniendo en mente el simplificar los cálculos computacionales, se considera un elemento de referencia ˆ en los ejes coordenados (ε,η) cuyos vértices están el ( 1, 1),(1, 1),(1,1) y ( 1,1) respectivamente, en el cual mediante una 4849 función afín será proyectado cualquier elemento rectangular e cuyos vértices (x (e) 1,y(e) 1 ),(x(e) 2,y(e) 2 ),(x(e) 3,y(e) 3 ) y (x(e) 4,y(e) 4 ) están tomados en sentido contrario al movimiento de las manecillas del reloj como se muetra en la figura Figura 1: mediante la transformación f(x, y) = T(ε, η)+b, quedando dicha transformación como x = x(e) 2 y = x(e) 2 enlacuallamatrizt estádadapor T = 2 x (e) 1 ε+ y(e) 2 y (e) 1 4 x (e) 1 ε+ y(e) 4 y (e) 1 ( x (e) 2 x(e) 1 2 x (e) 4 x(e) y (e) 2 y(e) 1 2 y (e) 4 y(e) 1 2 ) η (131) η (132) y el vector b= (b 1,b 2 ) es la posición del vector centroide del rectángulo e, 4950 también se tiene que la transformación inversa es ε = η = x b 1 y(e) ( x (e) 4 x(e) y(e) 1 2 ( x (e) y b 2 ) ( x b 1 y(e) 4 x(e) 1 2 x (e) 2 x(e) y(e) 1 2 x (e) 2 x(e) 1 2 y b 2 ) x b 1 y(e) ) 2 y(e) 1 2 x (e) 2 x(e) y(e) 4 y(e) 1 2. (133) Entonceslasφ (e) i quedandefinidasentérminosde ˆφ i como ˆφ 1 (ε,η) = 1 (1 ε)(1 η) (134) 4 ˆφ 2 (ε,η) = 1 4 (1+ε)(1 η) ˆφ 3 (ε,η) = 1 4 (1+ε)(1+η) ˆφ 4 (ε,η) = 1 4 (1 ε)(1+η) ylasfuncionesφ (e) i sonobtenidasporelconjuntoφ (e) i (x,y)=ˆφ i (ε,η)con(x,y)y (ε,η) relacionadas por la Ec.(131), entonces se tendrian las siguientes integrales ( [ K (e) φ (e) φ (e) i j + φ(e) φ (e) ] ) i j ij = a +cφ (e) i φ (e) j dxdy (135) e x x y y ([ ( )( ˆφ = a i ε + ˆφ i η ˆφj ε + ˆφ ) j η + ˆ ε x η x ε x η x ( )( ˆφ i ε + ˆφ i η ˆφj ε + ˆφ ) j η )]+cˆφ ε y η y ε y η y iˆφj J dεdη dondeelíndiceiyjvariade1a4.enestáúltimausamoslaregladelacadena ydxdy= J dεdηparaelcambiodevariableenlasintegrales,aquí J =dett, dondet estádadocomoenlaec.(132).pararesolver e f φ (e) j dxdyencada e segeneralasintegrales F (e) j = = dondeelíndiceiyj variade1a4. f φ (e) j dxdy (136) e f ˆφj J dεdη ˆ 5051 Para realizar el cálculo numérico de las integrales en el rectángulo de referencia ˆ=[ 1,1] [ 1,1],debemosconocer φi ε, φi η, ε x, ε y, η x y η y,entonces realizando las operaciones necesarias a la Ec.(134) obtenemos φ 1 ε = 1 4 (1 η) φ 1 η = 1 4 (1 ε) φ 2 ε = 1 4 (1 η) φ 2 η = 1 4 (1+ε) φ 3 ε = 1 4 (1+η) φ 3 η = 1 4 (1+ε) φ 4 ε = 1 4 (1+η) φ 4 η = 1 4 (1 ε) (137) ytambién ( ε x = ( η x = y (e) 4 y(e) 1 2detT y (e) 2 y(e) 1 2detT ) ) ( ε y = ( η y = x (e) 4 x(e) 1 2detT x (e) 2 x(e) 1 2detT ) ) (138) las cuales deberán de ser sustituidas en cada K (e) ij y F (e) j para calcular las integralesenelelemento e.estasintegralesseharánenelprogramausando cuadratura Gaussiana, permitiendo reducir el número de cálculos al mínimo pero manteniendo el balance entre precisión y número bajo de operaciones necesarias para realizar las integraciones. Suponiendo que fue dividido en E elementos, estos elementos generan N nodos en total, de los cuales N d son nodos desconocidos y N c son nodos conocidosconvalorγ j,entonceselalgoritmodeensambledelamatrizk yel vector F se puede esquematizar como: K i,j =(φ i,φ j ) i=1,2,...,e,j=1,2,...,e F j =(f,φ j ) j=1,2,...,e j=1,2,...,n d : b j =b j γ i K i,j i=1,2,...,e Así, se construye una matriz global en la cual están representados los nodos conocidos y los desconocidos, tomando sólo los nodos desconocidos de la matriz KformaremosunamatrizA,haciendolomismoalvectorF formamoselvector b, entonces la solución al problema será la resolución del sistema de ecuaciones lineales Ax= b, este sistema puede resolverse usando por ejemplo el método de Gradiente Conjugado. El vector x contendrá la solución buscada en los nodos desconocidosn d Método de Elemento Finito Usando Discretización de Triángulos Para resolver la Ec.(73), usando una discretización con triángulos, primero dividimoseldominio R 2 enn x nodoshorizontalesporn y nodosverticales, teniendoe=2(n x 1)(N y 1)subdominiosoelementostriángulares e tales 5152 que= E e=1 e y i j sisonadyacentes,conuntotalden =N x N y nodos. Dondelasfuncioneslinealesdefinidasporpedazosen e ennustrocasoserán polinomios de orden uno en cada variable separadamente y cuya restricción de φ i a e esφ (e) i (.Parasimplificarloscálculosenestaetapa,supondremosquela ) 1 0 matriza=a,entoncessetienequelaintegraldelladoizquierdode 0 1 la Ec.(174) queda escrita como ( ) a φi φ j +cφ i φ j dxdy= f φ j dxdy (139) donde K ij = = = E e=1 E e=1 ( a φi φ j +cφ i φ j ) dxdy (140) e e ( ) a φ (e) i φ (e) j +cφ (e) i φ (e) j dxdy ( a [ φ (e) φ (e) i j x + φ(e) i x y φ (e) ] ) j +cφ (e) i φ (e) j dxdy y yelladoderechocomo F j = = E e=1 f φ j dxdy (141) e f φ (e) j dxdy. Paracada e delasubmatrizdeintegrales(matrizdecargalocal) K ij = e ( tiene la estructura a [ φ (e) φ (e) i j x + φ(e) i x y k (e) 1,1 k (e) 1,2 k (e) 1,3 k (e) 2,1 k (e) 2,2 k (e) 2,3 k (e) 3,1 k (e) 3,2 k (e) 3,3 φ (e) ] ) j +cφ (e) i φ (e) j dxdy (142) y lacualdeberáserensambladaenlamatrizdecargaglobalquecorrespondaa lanumeracióndenodos localesdel elemento e conrespecto alanumeración globaldeloselementosen. De manera parecida, para cada e de se genera el vector de integrales (vector de carla gocal) F j = e f φ (e) j dxdy (143) 5253 con la estructura F (e) 1 F (e) 2 F (e) 3 el cual también deberá ser ensamblado en el vector de carga global que corresponda a la numeriación de nodos locales al elemento e con respecto a la numeración global de los elementos de. Montando los K (e) ij en la matriz K y los F (e) j en el vector F según la numeracióndenodosglobal,segeneraelsistemaku h =Fdondeu h seráelvector cuyos valores serán la solución aproximada a la Ec.(73) en los nodos interiores de.lamatrizkgeneradadeestaforma,tieneunapropiedadmuyimportante, esbandadayelanchodebandaesde7elementos,estoesmuyútilalmomento de soportar la matriz en memoria. y Paraimplementarnuméricamenteencada e lasintegrales ( [ φ (e) φ (e) i j + φ(e) φ (e) ] ) i j a +cφ (e) i φ (e) j dxdy (144) x x y y e e f φ (e) j dxdy teniendo en mente el simplificar los cálculos computacionales se considera a un elemento de referencia ˆ en los ejes coordenados (ε,η) cuyos vertices estanen(0,0),(1,0)y(0,1)yenelcualmedianteunmapeoafínseráproyectado caulquierelementotriangular e cuyosvertices(x (e) 1,y(e) 1 ),(x(e) 2,y(e) 2 ),(x(e) están tomados en el sentido contrario al movimiento de las manecillas del reloj como se muetra en la figura mediante la transformación f(ε,η) = T(ε,η)+b, quedando dicha transformación como yenlacuallamatrizt estádadapor 3,y(e) 3 ) x=x (e) 1 (1 ε η)+x(e) 2 ε+x(e) 3 η y=y (e) 1 (1 ε η)+y(e) 2 ε+y(e) 3 η (145) T = ( x (e) 2 x (e) 1 x (e) 3 x (e) 1 y (e) 2 y (e) 1 y (e) 3 y (e) 1 ) (146) dondebesunvectorconstante b= ( x (e) 1 y (e) 1 ) (147) 5354 Figura 2: también se tiene que la transformación inversa es donde ε = η = 1 [( )( y (e) 3 y (e) 1 2A e 1 [ 2A e ( y (e) 2 y (e) 1 x x (e) 1 )( A e = det ) x x (e) 1 ( )( x (e) 3 x (e) 1 ) 1 x (e) 1 y (e) 1 1 x (e) 2 y (e) 2 1 x (e) 3 y (e) 3 ( x (e) 2 x (e) 1 Entoceslasφ (e) i quedandefinidasentérminosdeˆφ i como y y (e) 1 )( )] y y (e) 1 )] (148). (149) ˆφ 1 (ε,η) = 1 ε η (150) ˆφ 2 (ε,η) = ε ˆφ 3 (ε,η) = η entoceslasfuncionesφ (e) i sonobtenidasporelconjuntoφ (e) i (x,y)=ˆφ i (ε,η)con (x, y) y(ε,η) relacionadas por la Ec.(145), entonces se tendrian las siguientes 5455 integrales k (e) ij = = e ˆ ( ([ a a [ φ (e) i x ( ˆφ i ε ε ( ˆφ i ε + ˆφ i ε y η φ (e) j + φ(e) i x y + ˆφ i x η η y η x )( ˆφj ε φ (e) j ] y )( ˆφj ε ε + ˆφ j y η +cφ (e) i φ (e) j ε + ˆφ j x η η y ) ) η + x ) )]+cˆφ iˆφj dxdy (151) J dεdη dondeelíndiceiyjvariade1a3.enestáúltimausamoslaregladelacadena ydxdy= J dεdηparaelcambiodevariableenlasintegrales,aquí J =dett, dondet estádadocomoenlaec.(146).pararesolver e f φ (e) j dxdyencada e segeneralasintegrales F (e) j = = f φ (e) j dxdy (152) e f ˆφj J dεdη ˆ dondeelindiceiyj varia1a3. Para realizar el cálculo numérico de las integrales en el triángulo de referencia ˆ, debemos conocer φi ε, φi η, ε x, ε y, η x y η y, entonces realizando las operaciones necesarias a las Ec.(150) obtenemos φ 1 ε = 1 φ 1 η = 1 φ 2 ε =1 φ 2 η =0 φ 3 ε =0 φ 3 η =1 (153) ytambién ( ) ε x = y (e) 3 y(e) 1 2A e ( η x = y (e) 2 y(e) 1 2A e ) ( ) ε y = x (e) 3 x(e) 1 2A e ( η y = x (e) 2 x(e) 1 ) (154) 2A e las cuales deberán de ser sustituidas en cada K (e) ij y F (e) j para calcular las integralesenelelemento e. Suponiendo que fue dividido en E elementos, estos elementos generan N nodos en total, de los cuales N d son nodos desconocidos y N c son nodos conocidosconvalorγ j,entonceselalgoritmodeensambledelamatrizk yel vector F se puede esquematizar como: K i,j =(φ i,φ j ) i=1,2,...,e,j=1,2,...,e F j =(f,φ j ) j=1,2,...,e j=1,2,...,n d : 5556 b j =b j γ i K i,j i=1,2,...,e Así, se construye una matriz global en la cual están representados los nodos conocidos y los desconocidos, tomando sólo los nodos desconocidos de la matriz KformaremosunamatrizA,haciendolomismoalvectorF formamoselvector b, entonces la solución al problema será la resolución del sistema de ecuaciones lineales Ax= b, este sistema puede resolverse usando por ejemplo el método de gradiente conjugado. El vector x contendrá la solución buscada en los nodos desconocidosn d. 5657 6. Métodos de Descomposición de Dominio(DDM) La solución numérica por los esquemas tradicionales de discretización tipo elemento finito y diferencias finitas generan una discretización del problema, la cual es usada para generar un sistema de ecuaciones algebraicos Au = b. Este sistema algebraico en general es de gran tamaño para problemas reales, al ser estos algoritmos secuenciales su implantación suele hacerse en equipos secuenciales y por ello no es posible resolver muchos problemas que involucren el uso de una gran cantidad de memoria, actualmente para tratar de subsanar dicha limitante, se usa equipo paralelo para soportar algoritmos secuenciales, haciendo ineficiente su implantación en dichos equipos. Los métodos de descomposición de dominio son un paradigma natural usado por la comunidad de modeladores. Los sistemas físicos son descompuestos en dos o más subdominios contiguos basados en consideraciones fenomenológicas. Estas descomposiciones basadas en dominios físicos son reflejadas en la ingeniería de software del código correspondiente. Los métodos de descomposición de dominio permiten tratar los problemas de tamaño considerable, empleando algoritmos paralelos en computadoras secuenciales y/o paralelas. Esto es posible ya que cualquier método de descomposición dedominiosebasaenlasuposicióndequedadoundominiocomputacional, estesepuedeparticionarensubdominios i, i=1,2,...,eentreloscualespuede o no existir traslape. Entonces el problema es reformulado en términos de cada subdominio(empleando algún método del tipo elemento finito) obteniendo una familia de subproblemas de tamaño reducido independientes en principio entre si,queestánacopladosatravésdelasoluciónenlainterfazdelossubdominios que es desconocida. Deestamanera,podemosclasificardemaneraburdaalosmétodosdedescomposición de dominio, como aquellos en que: existe traslape entre los subdominiosyenlosquenoexistetraslape.alaprimeraclaseperteneceelmétodo deschwarz(enelcualeltamañodeltraslapeesimportanteenlaconvergencia del método) y a los de la segunda clase pertenecen los métodos del tipo subestructuración(en el cual los subdominios sólo tienen en común los nodos de la frontera interior). La computación en paralelo es una técnica que nos permite distribuir una gran carga computacional entre muchos procesadores. Y es bien sabido que una de las mayores dificultades del procesamiento en paralelo es la coordinación de las actividades de los diferentes procesadores y el intercambio de información entre los mismos[12] mediante el paso de mensajes. Así, mediante los métodos de descomposición de dominio, la programación orientada a objetos y esquemas de paralelización que usan el paso de mensajes, es posible construir aplicaciones que coadyuven a la solución de problemas concomitantes en ciencia e ingeniería, ya que permiten utilizar todas las capacidades del cómputo en paralelo(supercomputadoras, clusters o grids), de esta forma es posible atacar una gran variedad de problemas que sin estas técnicas es imposible hacerlo de manera flexible y eficiente. Pero hay que notar que existen una amplia gama de problemas que nos 5758 interesan resolver que superan las capacidades de cómputo actuales, ya sea por el tiempo requerido para su solución, por el consumo excesivo de memoria o ambos. Lalistadelosmétodosdedescomposicióndedominioyeltipodeproblemas que pueden ser atacados por estos, es grande y está en constante evolución, ya que se trata de encontrar un equilibrio entre la complejidad del método (aunadaalapropiacomplejidaddelmodelo),laeficienciaenelconsumodelos recursos computacionales y la precisión esperada en la solución encontrada por los diversos métodos y las arquitecturas paralelas en la que se implante. A continuación describiremos algunos de estos métodos generales. En este capítulo se considerarán problemas con valor en la frontera(vbvp) de la forma Lu = f en (155) u = g en donde Lu= a u+cu (156) comouncasoparticulardeloperadorelípticodefinidoporlaec.(??)deorden dos Método de Schwarz El método fue desarrollado por Hermann Amandus Schwarz en 1869 (no como un método de descomposición de dominio), ya que en esos tiempos los matemáticos podían resolver problemas con geometrías sencillas de manera analítica, pero no tenían una idea clara de como poder resolver problemas que involucraran el traslape de esas geometrías sencillas. Como se conocía la solución para las geometrías sencillas por separado, la idea de Schwarz fue usar estas para conocer la solución en la geometría resultante al tener traslape, para más detalle ver[1]. Para describir el método, consideremos primero un dominio que está formadodedossubdominios 1 y 2 traslapados,esdecir 1 2,entonces = 1 2 ydenotemosaσ 1 = 1 2,Σ 2 = 2 1 y 1,2 = 1 2, como se muestra en la figura para dos dominios distintos: Figura3:Dominiosubdivididoendossubdominios 1 y 2. 5859 La forma original del método iterativo de Schwarz conocido como métodos alternantes de Schwarz, consiste en resolver sucesivamente los siguientes problemas. Sea u o una función de inicialización definida en, que se nulifica en, además hacemos u 0 1 =u 0 1 y u 0 2 =u 0 2. Parak 0definimos dos sucesiones u k+1 1 yu k+1 2 para resolver respectivamente y Lu k+1 1 =f en 1 u k+1 1 =u k 2 enσ 1 u k+1 1 =0 en 1 Lu k+1 2 =f en 2 u k+1 2 =u k+1 1 enσ 2 u k+1 2 =0 en 2 (157) (158) resolviendo los problemas secuencialmente en cada subdominio(por ejemplo con el método de elemento finito). Este método se conoce como Schwarz multiplicativo. ElmétodoalternantedeSchwarzdadoporlasEcs.(157)y(158)converge alasoluciónude(155)sisuponemosalgunasuavidadenlossubdominios 1 y 2,yaqueexistenconstantesC 1 yc 2 (0,1)talqueparatodok 0setiene u 1 u k+1 1 L ( 1 C k ) 1C2 k u u 0 L (Σ 1) (159) u 2 u k+1 2 C1 k+1 C2 k u u 0 L (Σ 2) L ( 2) lasconstantesc 1 yc 2 dereduccióndeerrordebendeestarbastantecercade 1 si la región de traslape 1,2 es delgada, la prueba de esta estimación puede encontrarseen[?]. Porotrolado,teniendoelconjuntou 0 1 =u0 1 yu 0 2 =u0 2,podemosgenerar dos pasos independientes uno de otro y Lu k+1 1 =f en 1 u k+1 1 =u k 2 enσ 1 u k+1 1 =0 en 1 Lu k+1 2 =f en 2 u k+1 2 =u k 1 enσ 2 u k+1 2 =0 en 2 (160) (161) resolviendo los problemas en paralelo de cada subdominio(por ejemplo con el método de elemento finito). Este método se conoce como Schwarz aditivo. La convergencia de este método en general requiere de algunas hipótesis adicionales, pero si converge, el número de iteraciones necesarias para converger será del doble que el método Schwarz multiplicativo. 5960 LageneralizacióndelmétododeSchwarzenelcasoenqueesparticionada en E > 2 subdominios traslapados puede describirse como sigue: DescomponiendoeldominioenEsubdominios e contraslapecomopor ejemplo, la descomposición siguiente Figura 4: Descomposición de en múltiples subdominios con traslape para el método de Schwarz. Entonces para resolver el problema por el método de Schwarz, primeramente esnecesariodefinirlossiguientessubespaciosdelespaciodesobolevh 1 (),en ellos estarán definidas las funciones usadas en el desarrollo del método: V i = { v i H 1 ( i ) v i =0 } V 0 i = H 1 0( i ) V i = { v H 1 0() v=0en\ i }. Denotaremos por I al operador identidad, por J i,i = 1,...E la inmersión de Vi sobre V (i.e. J i v = v para toda v Vi ), y por Ji T : H0( 1 i ) Vi al transpuestodeloperadorj i definidopor J T i F,v = F,J i v, F V,v Vi. (162) ydefinimos P i =J i P i :H 1 0( i ) H 1 0( i ). SeaL i :Vi 0 ( Vi 0 ) larestriccióndeloperadorlalsubespaciov 0 i,definido como L i w i,v i =a(w i,v i ) paratodaw i,v i Vi 0 ycomounoperadordeextensiónρ T i :V 0 i V i definido como ρ T iv i =ṽ i paratodav i V 0 i (163) 6061 yeltranspuestodeloperadorrestricciónρ i :(Vi ) (Vi 0) como ρ i G,v i = G,ρ T i v i,paratodag (V i ),v i Vi 0. (164) Deloanteriorsetieneque y L i =ρ i J T i LJ i ρ T i P i =J i P i :V V parai=1,...,e. (165) Entonces el método Multiplicativo de Schwarz queda como u k+ i E =(I Pi )u k+i 1 E +Ji ρ T i L 1 i ρ i J T i f (166) parai=1,2,...,e y ( E M u k+1 = I P i )u k + J i ρ T i L 1 i ρ i Ji T f (167) i=1 y la correspondiente ecuación de error como i=1 u u k+1 =(I P m )...(I P 1 )(u u k ). (168) ElelmétodoAditivodeSchwarzquedacomo ( ) E (u u u u k+1 k = I P ) i (169) i=1 y la correspondiente ecuación de error como Observaciones: u u k+1 =(I P m )...(I P 1 )(u u k ). (170) La precisión del método depende fuertemente del número de iteraciones realizadas en el proceso iterativo y converge a la precisión usadaenlasolucióndecadasubdominioenelmejordeloscasos. ElmétodoaditivodeSchwarzessecuencial,enelcasodelmétodo multiplicativo de Schwarz es paralelizable pero tiene una parte serial importante en el algoritmo y su convergencia no es la óptima en esta formulación, pero existen variantes del método que permiten remediaresto,paramásdetallesver[?],[4]y[5]. Hayquenotarqueporcadasubdominio(supóngasen)yencada iteración (supóngase I) se resuelve un problema para cada i, estosignificaquesiseusaelmétododeelementofinitopararesolver elproblemalocaldondeseusenenpromedior iteracionespararesolver el sistema lineal(no tomando en cuenta el costo invertido en generar las matrices), el total de iteraciones necesarias para resolver elproblemaeneldominioserár n I,resultandomuycostoso computacionalmente con respecto a otros métodos de descomposición de dominio. 6162 6.2. Método de Descomposición de Dominio de Subestructuración La solución numérica por los esquemas tradicionales de discretización tipo Elemento Finito y Diferencias Finitas generan una discretización del problema, la cual es usada para generar un sistema de ecuaciones algebraicos Au = b. Este sistema algebraico en general es de gran tamaño para problemas reales, al ser estos algoritmos secuenciales su implantación suele hacerse en equipos secuenciales y por ello no es posible resolver muchos problemas que involucren el uso de una gran cantidad de memoria, actualmente para tratar de subsanar dicha limitante, se usa equipo paralelo para soportar algoritmos secuenciales, haciendo ineficiente su implantación en dichos equipos. Los métodos de descomposición de dominio son un paradigma natural usado por la comunidad de modeladores. Los sistemas físicos son descompuestos en dos o más subdominios contiguos basados en consideraciones fenomenológicas. Estas descomposiciones basadas en dominios físicos son reflejadas en la ingeniería de software del código correspondiente. Los métodos de descomposición de dominio permiten tratar los problemas de tamaño considerable, empleando algoritmos paralelos en computadoras secuenciales y/o paralelas. Esto es posible ya que cualquier método de descomposición dedominiosebasaenlasuposicióndequedadoundominiocomputacional, estesepuedeparticionar-triangular-enesubdominios α, α=1,2,...,eentre los cuales no existe traslape. Entonces el problema es reformulado en términos de cada subdominio(empleando por ejemplo algún método tipo elemento finito) obteniendo una familia de subproblemas de tamaño reducido independientes en principioentresi,queestánacopladosatravésdelasoluciónenlainterfazdelos subdominios que es desconocida[2],[1],[5] y[7], uno de los primeros Métodos de Descomposición de Dominio sin traslape es el de Subestructuración-mejor conocidocomoschur-yelcualeslabaselosmétodosmáscomúnmenteusados como son FETI-DP, BDDC y el propio DVS-DDM. En esta sección y sin perdida de generalidad se considerarán problemas con valorenlafrontera(vbvp)delaforma donde Lu = f en (171) u = g en Lu= a u+cu (172) conaunamatrizpositivadefinida,simétricayc 0,comouncasoparticulardel operadorelípticodeorden2yparaejemplificartomaremosundominio R 2 con fronteras poligonales, es decir, es un conjunto abierto acotado y conexo talquesufrontera eslaunióndeunnúmerofinitodepolígonos. Simultiplicamosalaecuación a u+cu=f porv V =H 1 0 (), obtenemos v ( a u+cu ) =vf (173) 6263 aplicando el teorema de Green obtenemos ( ) v a u+cuv dx= vf dx. (174) Definiendo el operador bilineal ( ) a(u,v)= v a u+cuv dx (175) y la funcional lineal l(v)= f,v = vf dx (176) entonces podemos reescribir el problema en forma variacional, haciendo uso de laformabilineala(, )ylafuncionallineall( ). Donde las funciones lineales definidas por pedazos en e en nuestro caso serán polinomios de orden uno en cada variable separadamente y cuya restricción deφ i a e esφ (e) i (. Para simplificar ) los cálculos en esta etapa, supondremos que 1 0 lamatriza=a,entoncessetienequelaintegraldelladoizquierdo 0 1 delaec.(174)quedaescritacomo ( ) a φi φ j +cφ i φ j dxdy= f φ j dxdy (177) donde K ij = = = E e=1 E e=1 ( a φi φ j +cφ i φ j ) dxdy (178) ( ) a φ e (e) i φ (e) j +cφ (e) i φ (e) j dxdy e yelladoderechocomo ( a [ φ (e) φ (e) i j x F j = = + φ(e) i x y E e=1 φ (e) ] ) j +cφ (e) i φ (e) j dxdy y f φ j dxdy (179) e f φ (e) j dxdy. ConsideremoselproblemadadoporlaEc.(171)eneldominio,elcuales subdivididoenesubdominios i,i=1,2,...,esintraslape,tambiénconocida comomallagruesat H,esdecir i j = i j y = 63 E i, (180) i=164 y al conjunto Γ= E Γ i, siγ i = i \ (181) i=1 lo llamaremos la frontera interior del dominio, denotamos por H al diámetro H i = Diam( i ) de cada i que satisface Diam( i ) H para cada i = 1,2,...,E, además, cada subdominio i es descompuesto en una mallado finot h dek subdominiosmedianteunatriangulación e demodoqueestasea conforme,denotamosporhaldiámetroh i =Diam( e )decada e quesatisface Diam( e ) hparacadae=1,2,...,k decadai=1,2,...,e. Unejemplodeundominioysudescomposiciónensubdominios i ycada i asuvezdescompuestoen e subdominiossemuestraenlafigura: Figura 5: Dominio descompuesto en una partición gruesa de 3 3 y cada subdominio i enunaparticiónfinade7 5. Sinperdidadegeneralidadtomemosg=0en,notemosquesiemprees posibleponerelproblemadelaec.(171)comounoconcondicionesdefrontera Dirichlet que se nulifiquen mediante la adecuada manipulación del término del lado derecho de la ecuación. PrimeramenteseaD H 1 0 ()unespaciolinealdefuncionesdedimensión finitan,enelcualestédefinidounproductointeriordenotadoparacadau,v Dpor u v= u,v (182) Considerando la existencia de los subconjuntos linealmente independientes los cuales satisfacen B D,B I D I,B Γ D Γ B Γ D Γ,B ΓJ D Γ1,B ΓM D Γ2 (183) B=B I B Γ y B Γ =B ΓJ B ΓM (184) 6465 elespaciogeneradoporcadaunodelossubconjuntosb Γ y B Γ es D Γ,sinembargo distinguimoslapropiedaddequelosmiembrosdeb Γ tienensoportelocal. Definimos las bases { B I = w 1 I,...,w N I I },B ΓM = delasfuncionaleslinealesφ i en. { w 1 M,...,w N Γ M } { yb ΓJ = w 1 J,...,w N Γ J } (185) Entoncesdefiniendoparatodaδ=1,...,K,lamatrizdeN δ N δ ] A [ w δ i II I,wj I (186) quesóloestadefinidaencadasubespacio(subdominio δ ).Entonces,lamatriz virtuala II esdadaporlamatrizdiagonaldelaforma A II A 1 II A 2 II... A E II dondeelrestodelamatrizfueradeladiagonalenbloquesescero. De forma similar definimos A δ [ w i ] IΓ I,wα Γ, A δ [ w α ] ΓI Γ,wi I y (187) (188) A δ ΓΓ [ wα Γ,w α Γ ] (189) paratodaδ=1,...,e,obsérvesequecomo B Γ =B ΓJ B ΓM entonces A δ ΓΓ =[ wα Γ,w α Γ ]= [ w α Γ J,w α Γ J ] + [ w α ΓM,w α Γ M ] (190) ( ) T. también que A δ = A δ IΓ ΓI Entonces las matrices virtuales AΓI,A y A IΓ ΓΓ quedarán definidas como A 1 IΓ A 2 A IΓ IΓ (191). A E IΓ [ ] A ΓI A 1 A 2 A E (192) ΓI ΓI ΓI y A ΓΓ [ E i=1 A i ΓΓ ] (193) 6566 [ E ] donde i=1 Ai esconstruidasumandolasa i segúnelordendelosnodos ΓΓ ΓΓ globales versus los nodos locales. También consideremos al vector u (u 1,..,u E ) el cual puede ser escrito comou=(u I,u Γ )dondeu I =(u 1,..., NI )yu Γ =(u 1,..., NΓ ). Así, el sistema virtual quedando expresado como A 1 II A 1 ΓI... A E II A E ΓI A II u I +A IΓ u Γ = b I (194) A ΓI u I +A ΓΓ u Γ = b Γ u I1. u IE u I1. u IE + + A 1 IΓ. A E IΓ A ΓΓ u Γ1. u ΓE u Γ1. u ΓE = = b I1. b IE b Γ1. b ΓE omáscompactamentecomoau=b,notemosquelasmatricesa i ΓΓ,Ai ΓI,Ai y IΓ A i sonmatricesbandadas. II Siahoradespejamosu I delaprimeraecuacióndelsistemadadoporlaec. (194) obtenemos ) 1 ) u I = (A II (b I A IΓ u Γ si sustituimos u I en la segunda ecuación del sistema dado por la Ec. (194) entonces tenemos (A ΓΓ A ΓI (A II ) 1AIΓ )u Γ =b Γ A ΓI ( A II ) 1bI (195) enlacuallosnodosinterioresnofiguranenlaecuaciónytodoquedaenfunción delosnodosdelafronterainterioru Γ. AlamatrizformadaporA ΓΓ A ΓI (A II ) 1AIΓ seleconocecomoelcomplementodeschurglobalyseledenotacomo S=A ΓΓ A ΓI (A II ) 1AIΓ. (196) En nuestro caso, como estamos planteando todo en términos de subdominios i,coni=1,...,e,entonceslasmatricesa i ΓΓ,Ai ΓI,Ai IΓ yai quedandefinidas II de manera local, así que procedemos a definir el complemento de Schur local como ( ) 1A S i =A i ΓΓ Ai A i i (197) ΓI II IΓ 6667 adicionalmente definimos b i =b Γi A i ΓI ElsistemadadoporlaEc.(195)loescribimoscomo yquedadefinidodemaneravirtualapartirde [ E [ E ] S i ]u Γ = b i i=1 ( A i II) 1bIi. (198) Su Γ =b (199) i=1 (200) [ E [ donde i=1 i] S E y i=1 i] b podrían ser construida sumando las S i y b i respectivamente según el orden de los nodos globales versus los nodos locales. El sistema lineal virtual obtenido de la Ec.(199) se resuelve-dependiendo de si es o no simétrico- eficientemente usando el método de Gradiente Conjugado o algunavariantedegmres,paraellonoesnecesarioconstruirlamatrizscon lascontribucionesdecadas i correspondientesalsubdominioi,loquehacemos espasaracadasubdominioelvectoru Γi correspondientealai-ésimaiteración delmétodoiterativousadoparaqueencadasubdominioseevaluéũ i Γ =S i u Γi localmenteyconelresultadoseformaelvectorũ Γ = E i=1ũγ i ysecontinuécon los demás pasos del método. Esto es ideal para una implementación en paralelo del método de Gradiente Conjugado o variante de GMRES. Unavez resuelto el sistema de la Ec. (200) enel que hemos encontradola solución para los nodos de la frontera interior u Γ, entonces debemos resolver localmentelosu Ii correspondientesalosnodosinterioresparacadasubespacio i,paraestoempleamos u Ii = ( ) 1 ) A (b i II Ii A iiγ u Γi (201) paracadai=1,2,...,e,quedandoasíresueltoelproblemaau=btantoenlos nodosinterioresu Ii comoenlosdelafronterainterioru Γi correspondientesa cadasubespacio i. ( ) 1 Observación37 NotemosquenormalmentelasmatriceslocalesS i y A i II no se construyen, ya que estas serian matrices densas y su construcción es computacionalmentemuycostosa,ycomosólonosinteresaelproductosy Γ,omás [ E precisamente i=1 i] S y Γ,entoncessillamamosy Γi alvectorcorrespondiente al subdominio i, entonces tendremos ( ( ) 1A ) ũ i Γ = A i ΓΓ Ai A i i y ΓI II IΓ Γi. (202) 6768 Para evaluar eficientemente esta expresión, realizamos las siguientes operaciones equivalentes x1 = A i y ΓΓ Γ i (203) ( ( ) 1A ) x2 = A i A i i y ΓI II IΓ Γi ũ Γ i = x1 x2 la primera y tercera expresión no tienen ningún problema en su evaluación, para la segunda expresión tendremos que hacer x3=a i IΓ y Γ i (204) con este resultado intermedio deberíamos calcular ( 1x3 x4= A II) i (205) pero como no contamos con A i II obteniendo ( A i II) 1,entoncesmultiplicamoslaexpresiónpor ( 1x3 A i II x4=ai A i II II) (206) al simplificar, tenemos A i x4=x3. (207) II Esta última expresión puede ser resuelta usando Factorización LU, Gradiente Conjugado o alguna variante de GMRES (cada una de estas opciones tiene ventajas y desventajas computacionales que deben ser evaluadas al momento de implementar el código para un problema particular). Una vez obtenido x4, podremos calcular x2=a i x4 (208) ΓI así ũ Γ i =x1 x2 (209) completandolasecuenciadeoperacionesnecesariaparaobteners i y Γi. Observación38 Enelcasodelcálculode ( 1bIi b i =b Γi A i A i ΓI II) (210) algo análogo al comentario anterior deberá de hacerse, ya que nuevamente está ( 1,porellodeberemosdeusarelsiguienteprocedimientopara involucrado A II) i evaluar eficientemente esta expresión, realizando las operaciones equivalentes y1= ( A i II) 1bIi (211) 6869 multiplicandopora i II alaúltimaexpresión,obtenemos A i II y1=ai II ( A i II) 1bIi (212) simplificando, tenemos ( ) A i y1=b II Ii (213) donde esta última expresión puede ser resuelta usando Factorización LU, Gradiente Conjugado o alguna variante de GMRES, luego hacemos y para finalizar el cálculo, calculamos y2=a i y1 (214) ΓI b i =b Γi y2. (215) Observación 39 En la evaluación de ( ) 1 ) u Ii = A (b i II Ii A iiγ u Γi (216) ( 1,porellodeberemosdeusarelsiguiente esta nuevamente involucrado A II) i procedimiento para evaluar eficientemente esta expresión, realizando las operaciones equivalentes x4 = b Ii A i u IΓ Γ i (217) ( 1x4 u Ii = A II) i multiplicandopora i II alaúltimaexpresión,obtenemos simplificando, tenemos A i II u I i =A i II ( A i II) 1x4 (218) A i II u I i =x4 (219) esta última expresión puede ser resuelta usando Factorización LU o Gradiente Conjugado. Comoseindicoenlasúltimasobservaciones,pararesolverelsistemaA i II x= b podemos usar Factorización LU, Gradiente Conjugado, alguna variante de GMRES o cualquier otro método para resolver sistemas lineales, pero deberá de usarse aquel que proporcione la mayor velocidad en el cálculo o que consuma la menor cantidad de memoria (ambas condicionantes son mutuamente excluyentes), por ello la decisión de que método usar deberá de tomarse al momento de tener que resolver un problema particular en un equipo dado y básicamenteelcondicionanteeseltamañodellamatriza i II. 6970 Para usar el método de Factorización LU, se deberá primeramente de factorizar la matriz bandada A i en una matriz LU, la cual es bandada pero II incrementaeltamañodelabandaamásdeldoble,peroestaoperaciónsólose deberá de realizar una vez por cada subdominio, y para solucionar los diversos sistemaslinealesa i x=bsóloseránecesarioevaluarlossistemas II Ly = b (220) Ux = y en donde y es un vector auxiliar. Esto proporciona una manera muy eficiente de evaluar el sistema lineal pero el consumo en memoria para un problema particular puede ser excesivo. Porello,sielproblemainvolucraunagrancantidaddenodosinterioresyel equipoenelqueseimplantarálaejecucióndelprogramatieneunacantidadde memoria muy limitada, es recomendable usar el método de Gradiente Conjugado o alguna variante de GMRES, este consume una cantidad de memoria adicional muy pequeña y el tiempo de ejecución se optimiza versus la Factorización LU. De esta forma, es posible adaptar el código para tomar en cuenta la implementaciónde este en unequipo de cómputo enparticular y poder sacar el máximo provecho al método de Subestructuración en la resolución de problemas elípticos de gran envergadura. Enloquerestadelpresentetrabajo,seasumequeelmétodoempleadopara resolvera i ( x=bensusrespectivasvariantesnecesariasparaevitarelcálculo II 1 de A II) i es el método de Gradiente Conjugado, logrando así el máximo desempeño en velocidad en tiempo de ejecución. El número de condicionamiento del complemento de Schur sin precondicionamiento puede ser estimado, para ello: Definición40 Introducimosunanorma-L 2 equivalentesobreγmediante u Γ 2 E Γ = u Γ 2 L 2 ( i ). (221) i=1 Teorema41 Sea u Γ la traza de funciones de elemento finito en V h sobre Γ, asumimos que los coeficientes de la ecuación diferencial parcial ρ i = 1, i = 1,2,...,E, y que la malla fina T h y la malla gruesa T H sea cuasi-uniforme. Entonces existen dos constantes positivas c y C, independientes de h y H, tal que ch uγ 2 Γ s(u Γ,u Γ ) Ch 1 uγ 2 (222) Γ deestemodo κ=cond(s) C Hh. (223) 7071 Por analogía al método de subestructuración desarrollado anteriormente, dado un sistema lineal Mx = f que proviene de la discretización de algún método tipo Diferencias Finitas, Elemento Finito o Volumen Finito, siempre es posible reacomodarlo como ( )( ) ( ) A B u a = (224) C D v b ( ) ( ) ( ) A B u a con M =, x= y f =, enla cual la matriz A sea C D v b invertible,entonces ( D CA 1 B ) v=b CA 1 a (225) ydonde u=a 1( a Bv ). (226) Así,hemostransformadoelsistemaMx=f,enotroequivalente peroconunmenornúmerodegradosdelibertad,donde Nv=g (227) N = ( D CA 1 B ), g=b CA 1 a. (228) a esta descomposición matricial se le conoce como la descomposición de Schur. Así, para solucionar el sistema lineal N v = g, seleccionaremos el método adecuado acorde a las propiedades de la matriz N como se vio enel presente capítulo, siendo los más comunes el método CGM-para matrices simétricas- y variantes del método GMRES-para matrices no simétricas-, entre otros. 7172 7. El Cómputo en Paralelo Los sistemas de cómputo con procesamiento en paralelo surgen de la necesidad de resolver problemas complejos en un tiempo razonable, utilizando las ventajas de memoria, velocidad de los procesadores, formas de interconexión deestosydistribucióndelatarea,alosqueensuconjuntodenominamosarquitectura en paralelo. Entenderemos por una arquitectura en paralelo a un conjunto de procesadores interconectados capaces de cooperar en la solución de un problema. Así, para resolver un problema en particular, se usa una o combinación de múltiples arquitecturas(topologías), ya que cada una ofrece ventajas y desventajas que tienen que ser sopesadas antes de implementar la solución del problema en una arquitectura en particular. También es necesario conocer los problemas a los que se enfrenta un desarrollador de programas que se desean correr en paralelo, como son: el partir eficientemente un problema en múltiples tareas y como distribuir estas según la arquitectura en particular con que se trabaje Arquitecturas de Software y Hardware En esta sección se explican en detalle las dos clasificaciones de computadoras más conocidas en la actualidad. La primera clasificación, es la clasificación clásicadeflynnendóndesetienenencuentasistemasconunoovariosprocesadores, la segunda clasificación es moderna en la que sólo tienen en cuenta los sistemas con más de un procesador. El objetivo de esta sección es presentar de una forma clara los tipos de clasificación que existen en la actualidad desde el punto de vista de distintos autores, así como cuáles son las ventajas e inconvenientes que cada uno ostenta, ya que es común que al resolver un problema particular se usen una o más arquitecturas de hardware interconectadas generalmente por red Clasificación de Flynn Clasificación clásica de arquitecturas de computadoras que hace alusión a sistemas con uno o varios procesadores, Flynn la publicó por primera vez en 1966yporsegundavezen1970. Estataxonomíasebasaenelflujoquesiguenlosdatosdentrodelamáquina y de las instrucciones sobre esos datos. Se define como flujo de instrucciones al conjunto de instrucciones secuenciales que son ejecutadas por un único procesadorycomoflujodedatosalflujosecuencialdedatosrequeridosporelflujo de instrucciones. Con estas consideraciones, Flynn clasifica los sistemas en cuatro categorías: Single Instruction stream, Single Data stream (SISD) Los sistemas de este tipo se caracterizan por tener un único flujo de instrucciones sobre un único flujo de datos, es decir, se ejecuta una instrucción detrás de otra. Este es el concepto de arquitectura serie de Von Neumann donde, en cualquier 7273 momento, sólo se ejecuta una única instrucción, un ejemplo de estos sistemas son las máquinas secuenciales convencionales. Single Instruction stream, Multiple Data stream (SIMD) Estos sistemas tienen un único flujo de instrucciones que operan sobre múltiples flujos de datos. Ejemplos de estos sistemas los tenemos en las máquinas vectoriales con hardware escalar y vectorial. El procesamiento es síncrono, la ejecución de las instrucciones sigue siendo secuencial como en el caso anterior, todos los elementos realizan una misma instrucción pero sobre una gran cantidad de datos. Por este motivo existirá concurrencia de operación, es decir, esta clasificación es el origen de la máquina paralela. Elfuncionamientodeestetipodesistemaseselsiguiente.Launidaddecontrol manda una misma instrucción a todas las unidades de proceso(alus). Las unidades de proceso operan sobre datos diferentes pero con la misma instrucción recibida. Existen dos alternativas distintas que aparecen después de realizarse esta clasificación: Arquitectura Vectorial con segmentación, una CPU única particionada en unidades funcionales independientes trabajando sobre flujos de datos concretos. Arquitectura Matricial(matriz de procesadores), varias ALUs idénticas a las que el procesador da instrucciones, asigna una única instrucción pero trabajando sobre diferentes partes del programa. Multiple Instruction stream, Single Data stream (MISD) Sistemas con múltiples instrucciones que operan sobre un único flujo de datos. Este tipo de sistemas no ha tenido implementación hasta hace poco tiempo. Los sistemas MISD se contemplan de dos maneras distintas: Varias instrucciones operando simultáneamente sobre un único dato. Varias instrucciones operando sobre un dato que se va convirtiendo en un resultado que será la entrada para la siguiente etapa. Se trabaja de forma segmentada, todas las unidades de proceso pueden trabajar de forma concurrente. Multiple Instruction stream, Multiple Data stream(mimd) Sistemas con un flujo de múltiples instrucciones que operan sobre múltiples datos. Estos sistemas empezaron a utilizarse antes de la década de los 80s. Son sistemas con memoria compartida que permiten ejecutar varios procesos simultáneamente (sistema multiprocesador). Cuando las unidades de proceso reciben datos de una memoria no compartida estos sistemas reciben el nombre de MULTIPLE SISD (MSISD). En 7374 arquitecturas con varias unidades de control(misd Y MIMD), existe otro nivel superior con una unidad de control que se encarga de controlar todas las unidades de control del sistema (ejemplo de estos sistemas son las máquinas paralelas actuales) Categorías de Computadoras Paralelas Clasificación moderna que hace alusión única y exclusivamente a los sistemas que tienen más de un procesador(i.e máquinas paralelas). Existen dos tipos de sistemas teniendo en cuenta su acoplamiento: Los sistemas fuertemente acoplados son aquellos en los que los procesadores dependen unos de otros. Los sistemas débilmente acoplados son aquellos en los que existe poca interacción entre los diferentes procesadores que forman el sistema. Atendiendo a esta y a otras características, la clasificación moderna divide a los sistemas en dos tipos: Sistemas multiprocesador(fuertemente acoplados) y sistemas multicomputadoras(débilmente acoplados). Multiprocesadores o Equipo Paralelo de Memoria Compartida Un multiprocesador puede verse como una computadora paralela compuesta por varios procesadores interconectados que comparten un mismo sistema de memoria. Los sistemas multiprocesadores son arquitecturas MIMD con memoria compartida. Tienen un único espacio de direcciones para todos los procesadores y losmecanismosdecomunicaciónsebasanenelpasodemensajesdesdeelpunto de vista del programador. Dado que los multiprocesadores comparten diferentes módulos de memoria, pudiendo acceder a un mismo módulo varios procesadores, a los multiprocesadores también se les llama sistemas de memoria compartida. Figura 6: Arquitectura de una computadora paralela con memoria compartida Para hacer uso de la memoria compartida por más de un procesador, se requiere hacer uso de técnicas de semáforos que mantienen la integridad de la memoria; esta arquitectura no puede crecer mucho en el número de procesadores interconectados por la saturación rápida del bus o del medio de interconexión. Dependiendo de la forma en que los procesadores comparten la memoria, se clasifican en sistemas multiprocesador UMA, NUMA, COMA y Pipeline. 7475 Uniform Memory Access(UMA) Sistema multiprocesador con acceso uniforme a memoria. La memoria física es uniformemente compartida por todos los procesadores, esto quiere decir que todos los procesadores tienen el mismo tiempodeaccesoatodaslaspalabrasdelamemoria.cadaprocesadortienesu propia caché privada y también se comparten los periféricos. Los multiprocesadores son sistemas fuertemente acoplados(tightly-coupled), dadoelaltogradodecomparticióndelosrecursos(hardwareosoftware)yelalto nivel de interacción entre procesadores, lo que hace que un procesador dependa deloquehaceotro. El sistema de interconexión debe ser rápido y puede ser de uno de los siguientes tipos: bus común, red crossbar y red multietapa. Este modelo es conveniente para aplicaciones de propósito general y de tiempo compartido por varios usuarios, existen dos categorías de sistemas UMA. Sistema Simétrico Cuando todos los procesadores tienen el mismo tiempo de acceso a todos los componentes del sistema (incluidos los periféricos), reciben el nombre de sistemas multiprocesador simétrico. Los procesadores tienen el mismo dominio(prioridad) sobre los periféricos y cada procesador tiene la misma capacidad para procesar. Sistema Asimétrico Los sistemas multiprocesador asimétrico, son sistemas con procesadores maestros y procesadores esclavos, en donde sólo los primeros pueden ejecutar aplicaciones y dónde en tiempo de acceso para diferentes procesadores no es el mismo. Los procesadores esclavos (attached) ejecutan código usuario bajo la supervisión del maestro, por lo tanto cuando una aplicación es ejecutada en un procesador maestro dispondrá de una cierta prioridad. Non Uniform Memory Access(NUMA) Un sistema multiprocesador NUMAesunsistemadememoriacompartidadondeeltiempodeaccesovaría según donde se encuentre localizado el acceso. El acceso a memoria, por tanto, no es uniforme para diferentes procesadores, existen memorias locales asociadas a cada procesador y estos pueden acceder a datos de su memoria local de una manera más rápida que a las memorias de otros procesadores, debido a que primero debe aceptarse dicho acceso por el procesador del que depende el módulo de memoria local. Todas las memorias locales conforman la memoria global compartida y físicamente distribuida y accesible por todos los procesadores. Cache Only Memory Access (COMA) Los sistemas COMA son un casoespecialdelossistemasnuma.estetipodesistemasnohatenidomucha trascendencia, al igual que los sistemas SIMD. 7576 Las memorias distribuidas son memorias cachés, por este motivo es un sistema muy restringido en cuanto a la capacidad de memoria global. No hay jerarquía de memoria en cada módulo procesador. Todas las cachés forman un mismo espacio global de direcciones. El acceso a las cachés remotas se realiza a través de los directorios distribuidos de las cachés. Dependiendo de la red de interconexión utilizada, se pueden utilizar jerarquías en los directorios para ayudar a la localización de copias de bloques de caché. Procesador Vectorial Pipeline En la actualidad es común encontrar en un solo procesador los denominados Pipeline o Procesador Vectorial Pipeline del tipo MISD. En estos procesadores los vectores fluyen a través de las unidades aritméticas Pipeline. Las unidades constan de una cascada de etapas de procesamiento compuestas de circuitos que efectúan operaciones aritméticas o lógicas sobre el flujo de datos que pasan a través de ellas, las etapas están separadas por registros de alta velocidad usados para guardar resultados intermedios. Así la información que fluyeentrelasetapasadyacentesestábajoelcontroldeunrelojqueseaplicaa todos los registros simultáneamente. Multicomputadoras o Equipo Paralelo de Memoria Distribuida Los sistemas multicomputadoras se pueden ver como una computadora paralela en el cual cada procesador tiene su propia memoria local. En estos sistemas la memoria se encuentra distribuida y no compartida como en los sistemas multiprocesador. Los procesadores se comunican a través de paso de mensajes, ya queéstossólotienenaccesodirectoasumemorialocalynoalasmemoriasdel resto de los procesadores. Figura 7: Arquitectura de una computadora paralela con memoria distribuida La transferencia de los datosse realizaatravés de laredde interconexión que conecta un subconjunto de procesadores con otro subconjunto. La transferencia de unos procesadores a otros se realiza por múltiples transferencias entre procesadores conectados dependiendo del establecimiento de dicha red. 7677 Dado que la memoria está distribuida entre los diferentes elementos de proceso, estos sistemas reciben el nombre de distribuidos. Por otra parte, estos sistemas son débilmente acoplados, ya que los módulos funcionan de forma casi independiente unos de otros. Este tipo de memoria distribuida es de acceso lento porserpeticionesatravésde lared,peroesunaformamuy efectivadetener accesoaungranvolumendememoria. Equipo Paralelo de Memoria Compartida-Distribuida La tendencia actual en las máquinas paralelas es de aprovechar las facilidades de programación que ofrecen los ambientes de memoria compartida y la escalabilidad de las ambientes de memoria distribuida. En este modelo se conectan entre si módulos de multiprocesadores, pero se mantiene la visión global de la memoria apesardequeesdistribuida. Clusters El desarrollo de sistemas operativos y compiladores del dominio público (Linux y software GNU), estándares para el pase de mensajes (MPI), conexión universal a periféricos(pci), etc. han hecho posible tomar ventaja de los económicos recursos computacionales de producción masiva (CPU, discos, redes). La principal desventaja que presenta a los proveedores de multicomputadorasesquedebensatisfacerunaampliagamadeusuarios,esdecir,debenser generales. Esto aumenta los costos de diseños y producción de equipos, así como los costos de desarrollo de software que va con ellos: sistema operativo, compiladores y aplicaciones. Todos estos costos deben ser añadidos cuando se hace una venta. Por supuesto alguien que sólo necesita procesadores y un mecanismo depasedemensajesnodeberíapagarportodosestosañadidosquenuncausará. Estos usuarios son los que están impulsando el uso de clusters principalmente de computadoras personales(pc), cuya arquitectura se muestra a continuación: Figura 8: Arquitectura de un cluster Los cluster se pueden clasificar en dos tipos según sus características físicas: Cluster homogéneo si todos los procesadores y/o nodos participantes en el equipo paralelo son iguales en capacidad de cómputo (en la cual es permitido variar la cantidad de memoria o disco duro en cada procesador). 7778 Clusterheterogéneoesaquelenquealmenosunodelosprocesadoresy/o nodos participantes en el equipo paralelo son de distinta capacidad de cómputo. Los cluster pueden formarse de diversos equipos; los más comunes son los de computadoras personales, pero es creciente el uso de computadoras multiprocesador de más de un procesador de memoria compartida interconectados por red con los demás nodos del mismo tipo, incluso el uso de computadoras multiprocesador de procesadores vectoriales Pipeline. Los cluster armados con la configuración anterior tienen grandes ventajas para procesamiento paralelo: Larecienteexplosiónenredesimplicaquelamayoríadeloscomponentes necesarios para construir un cluster son vendidos en altos volúmenes y por lo tanto son económicos. Ahorros adicionales se pueden obtener debido a que sólo se necesitará una tarjeta de vídeo, unmonitoryuntecladoporcluster.elmercadodelosmultiprocesadores es más reducido y más costoso. Remplazar un componente defectuoso en un cluster es relativamente trivial comparado con hacerlo en un multiprocesador, permitiendo una mayor disponibilidad de clusters cuidadosamente diseñados. Desventajas del uso de clusters de computadoras personales para procesamiento paralelo: Con raras excepciones, los equipos de redes generales producidos masivamente no están diseñados para procesamiento paralelo y típicamente su latencia es alta y los anchos de banda pequeños comparados con multiprocesadores. Dado que los clusters explotan tecnología que sea económica, los enlaces en el sistema no son veloces implicando que la comunicación entre componentes debe pasar por un proceso de protocolos de negociación lentos, incrementando seriamentelalatencia.enmuchosyenelmejordeloscasos(debido a costos) se recurre a una red tipo Fast Ethernet restringimiento la escalabilidad del cluster. Hay poco soporte de software para manejar un cluster como un sistema integrado. Los procesadores no son tan eficientes como los procesadores usados en los multiprocesadores para manejar múltiples usuarios y/o procesos. Esto hace que el rendimiento de los clusters se degrade con relativamente pocos usuarios y/o procesos. Muchas aplicaciones importantes disponibles en multiprocesadores y optimizadas para ciertas arquitecturas, no lo están en clusters. Sin lugar a duda los clusters presentan una alternativa importante para varios problemas particulares, no sólo por su economía, si no también porque 7879 pueden ser diseñados y ajustados para ciertas aplicaciones. Las aplicaciones que pueden sacar provecho de clusters son en donde el grado de comunicación entre procesosesdebajoamedio. Tipos de Cluster Básicamente existen tres tipo de clusters, cada uno de ellos ofrece ventajas y desventajas, el tipo más adecuado parael cómputo científicoes del de altorendimiento, pero existen aplicaciones científicas que pueden usar más de un tipo al mismo tiempo. Alta-disponibilidad (Fail-over o High-Availability): este tipo de cluster esta diseñado para mantener uno o varios servicios disponibles inclusoacostaderendimiento,yaquesufunciónprincipalesqueel servicio jamás tenga interrupciones como por ejemplo un servicio de bases de datos. Alto-rendimiento(HPC o High Performance Computing): este tipo de cluster está diseñado para obtener el máximo rendimiento de la aplicación utilizada incluso a costa de la disponibilidad del sistema, es decir el cluster puede sufrir caídas, este tipo de configuración está orientada a procesos que requieran mucha capacidad de cálculo. Balanceo de Carga(Load-balancing): este tipo de cluster esta diseñado para balancear la carga de trabajo entre varios servidores, lo que permite tener, por ejemplo, un servicio de cálculo intensivo multiusuarios que detecte tiempos muertos del proceso de un usuario para ejecutar en dichos tiempos procesos de otros usuarios. Grids Son cúmulos (grupo de clusters) de arquitecturas en paralelo interconectados por red, los cuales distribuyen tareas entre los clusters que lo forman, estos pueden ser homogéneos o heterogéneos en cuanto a los nodos componentes del cúmulo. Este tipo de arquitecturas trata de distribuir cargas de trabajo acorde a las características internas de cada cluster y las necesidades propiasdecadaproblema,estosehaceadosniveles,unaenlapartedeprogramaciónconjuntamenteconelbalancedecargasyotraenlapartedehardware que tiene que ver con las características de cada arquitectura que conforman al cúmulo Métricas de Desempeño Las métricas de desempeño del procesamiento de alguna tarea en paralelo es un factor importante para medir la eficiencia y consumo de recursos al resolver una tarea con un número determinado de procesadores y recursos relacionados de la interconexión de éstos. Entre las métricas para medir desempeño en las cuales como premisa se mantiene fijo el tamaño del problema, destacan las siguientes: Factor de aceleración, eficiencia y fracción serial. Cada una de ellas mide algo en particular 7980 y sólo la combinación de estas dan un panorama general del desempeño del procesamiento en paralelo de un problema en particular en una arquitectura determinada al ser comparada con otras. Factor de Aceleración(o Speed-Up) Se define como el cociente del tiempo que se tarda en completar el cómputo de la tarea usando un sólo procesador entre el tiempo que necesita para realizarlo en p procesadores trabajando en paralelo s= T(1) (229) T(p) en ambos casos se asume que se usará el mejor algoritmo tanto para un solo procesador como para p procesadores. Estamétricaenelcasoidealdeberíadeaumentardeformalinealalaumento del número de procesadores. Eficiencia Sedefinecomoelcocientedeltiempoquesetardaencompletarel cómputo de la tarea usando un solo procesador entre el número de procesadores multiplicado por el tiempo que necesita para realizarlo en p procesadores trabajando en paralelo e= T(1) pt(p) = s p. (230) Este valor será cercano a la unidad cuando el hardware se esté usando de manera eficiente, en caso contrario el hardware será desaprovechado. Fracción serial Sedefinecomoelcocientedeltiempoquesetardaencompletar el cómputo de la parte secuencial de una tarea entre el tiempo que se tarda el completar el cómputo de la tarea usando un solo procesador f= T s T(1) (231) perousandolaleydeamdahl T(p)=T s + T p p y rescribiéndola en términos de factor de aceleración, obtenemos la forma operativa del cálculo de la fracción serial que adquiere la forma siguiente f = 1 s 1 p 1 1 p. (232) Esta métrica permite ver las inconsistencias en el balance de cargas, ya que suvalordebieradetenderaceroenelcasoideal,porellounincrementoenel valor de f es un aviso de granularidad fina con la correspondiente sobrecarga en los procesos de comunicación. 8081 7.3. Cómputo Paralelo para Sistemas Continuos Como se mostró en los capítulos anteriores, la solución de los sistemas continuos usando ecuaciones diferenciales parciales genera un alto consumo de memoria e involucran un amplio tiempo de procesamiento; por ello nos interesa trabajar en computadoras que nos puedan satisfacer estas demandas. Actualmente, en muchos centros de cómputo es una práctica común usar directivas de compilación en equipos paralelos sobre programas escritos de forma secuencial, con la esperanza que sean puestos por el compilador como programas paralelos. Esto en la gran mayoría de los casos genera códigos poco eficientes, peseaquecorrenenequiposparalelosypuedenusartodalamemoriacompartida de dichos equipos, el algoritmo ejecutado continua siendo secuencial en la gran mayoría del código. Si la arquitectura paralela donde se implemente el programa es UMA de acceso simétrico, los datos serán accesados a una velocidad de memoria constante. Encasocontrario,alaccederaunconjuntodedatosescomúnqueunapartede estosseanlocalesaunprocesador(conunaccesodelordendenanosegundos), peroelrestodelosdatosdeberándeseraccesadosmediantered(conaccesodel orden de mili segundos), siendo esto muy costoso en tiempo de procesamiento. Por ello, si usamos métodos de descomposición de dominio es posible hacer que el sistema algebraico asociado pueda distribuirse en la memoria local de múltiples computadoras y que para encontrar la solución al problema se requiera poca comunicación entre los procesadores. Por lo anterior, si se cuenta con computadoras con memoria compartida o que tengan interconexión por bus, salvo en casos particulares no será posible explotar éstas características eficientemente. Pero en la medida en que se adecuen los programas para usar bibliotecas y compiladores acordes a las características del equipo disponible (algunos de ellos sólo existen de manera comercial) la eficiencia aumentará de manera importante. La alternativa más adecuada(en costo y flexibilidad), es trabajar con computadoras de escritorio interconectadas por red que pueden usarse de manera cooperativa para resolver nuestro problema. Los gastos en la interconexión de losequipossonmínimos(sóloelswitchyunatarjetaderedporequipoycables parasuconexión). Porellolosclustersylosgridssonenprincipiounabuena opción para la resolución de este tipo de problemas. Esquema de Paralelización Maestro-Esclavo La implementación de los métodos de descomposición de dominio que se trabajarán será mediante el esquema Maestro-Esclavo(Farmer) en el lenguaje de programación C++ bajo la interfaz de paso de mensajes MPI trabajando en un cluster Linux Debian. Donde tomando en cuenta la implementación en estrella del cluster, el modelo de paralelismo de MPI y las necesidades propias de comunicación del programa, el nodo maestro tendrá comunicación sólo con cada nodo esclavo y no existirá comunicación entre los nodos esclavos, esto reducirá las comunicaciones y optimizará el paso de mensajes. El esquema de paralelización maestro-esclavo, permite sincronizar por parte 8182 del nodo maestro las tareas que se realizan en paralelo usando varios nodos esclavos, éste modelo puede ser explotado de manera eficiente si existe poca comunicación entre el maestro y el esclavo y los tiempos consumidos en realizar las tareas asignadas son mayores que los períodos involucrados en las comunicaciones para la asignación de dichas tareas. De esta manera se garantiza que la mayoría de los procesadores estarán trabajando de manera continua y existirán pocos tiempos muertos. Figura 9: Esquema del maestro-esclavo Unfactorlimitanteenesteesquemaesqueelnodomaestrodeberádeatendertodaslas peticioneshechas porcada unode losnodosesclavos, esto toma especial relevancia cuando todos o casi todos los nodos esclavos compiten por ser atendidos por el nodo maestro. Se recomienda implementar este esquema en un cluster heterogéneo en donde el nodo maestro sea más poderoso computacionalmente que los nodos esclavos. Siaésteesquemaseleagregaunareddealtavelocidadydebajalatencia,sele permitirá operar al cluster en las mejores condiciones posibles, pero este esquema se verá degradado al aumentar el número de nodos esclavos inexorablemente. Perohayquesercuidadososencuantoalnúmerodenodosesclavosquese usan en la implementación en tiempo de ejecución versus el rendimiento general del sistema al aumentar estos, algunas observaciones posibles son: El esquema maestro-esclavo programado en C++ y usando MPI lanza P procesos (uno para el nodo maestro y P 1 para los nodos esclavos), estos en principio corren en un solo procesador pero pueden ser lanzados en múltiples procesadores usando una directiva de ejecución, de esta manera es posible que en una sola maquina se programe, depureysea puestoapunto el códigousandomallas pequeñas(delordendecientosdenodos)ycuandoestelistopuede mandarse a producción en un cluster. El esquema maestro-esclavo no es eficiente si sólo se usan dos procesadores(uno para el nodo maestro y otro para el nodo esclavo), yaqueelnodomaestroengeneralnorealizaloscálculospesadosy su principal función será la de distribuir tareas; los cálculos serán delegados al nodo esclavo. En el caso que nos interesa implementar, el método de descomposición de dominio adolece de este problema. 8283 Paso de Mensajes Usando MPI Para poder intercomunicar al nodo maestroconcadaunodelosnodosesclavosseusalainterfazdepasodemensajes (MPI), una biblioteca de comunicación para procesamiento en paralelo. MPI ha sido desarrollado como un estándar para el paso de mensajes y operaciones relacionadas. Este enfoque es adoptado por usuarios e implementadores de bibliotecas, en la cual se proveen a los programas de procesamiento en paralelo de portabilidad y herramientas necesarias para desarrollar aplicaciones que puedan usar el cómputo paralelo de alto desempeño. El modelo de paso de mensajes posibilita a un conjunto de procesos que tienen solo memoria local la comunicación con otros procesos(usando Bus o red) mediante el envío y recepción de mensajes. Por definición el paso de mensajes posibilita transferir datos de la memoria local de un proceso a la memoria local de cualquier otro proceso que lo requiera. En el modelo de paso de mensajes para equipos paralelos, los procesos se ejecutan en paralelo, teniendo direcciones de memoria separada para cada proceso, la comunicación ocurre cuando una porción de la dirección de memoria de unprocesoescopiadamedianteelenvíodeunmensajedentrodeotroproceso en la memoria local mediante la recepción del mismo. Las operaciones de envío y recepción de mensajes es cooperativa y ocurre sólo cuando el primer proceso ejecuta una operación de envío y el segundo proceso ejecuta una operación de recepción, los argumentos base de estas funciones son: Paraelqueenvía,ladireccióndelosdatosatransmitiryelproceso destino al cual los datos se enviarán. Parael que recibe, debedetenerladirecciónde memoriadonde se pondrán los datos recibidos, junto con la dirección del proceso del que los envío. Es decir: Send(dir, lg, td, dest, etiq, com) {dir, lg, td} describe cuántas ocurrencias lg de elementos del tipo de dato td se transmitirán empezando en la dirección de memoria dir. {des, etiq, com} describe el identificador etq de destino des asociado con la comunicación com. Recv(dir, mlg, td, fuent, etiq, com, st) {dir, lg, td} describe cuántas ocurrencias lg de elementos del tipo de dato td se transmitirán empezando en la dirección de memoria dir. {fuent, etiq, com, est} describe el identificador etq de la fuente fuent asociadoconlacomunicacióncom yelestadost. Elconjuntobásicodedirectivas(ennuestrocasosóloseusanestas)enC++ dempison: 8384 MPI::Init Inicializa al MPI MPI::COMM_WORLD.Get_size Busca el número de procesos existentes MPI::COMM_WORLD.Get_rank Busca el identificador del proceso MPI::COMM_WORLD.Send Envía un mensaje MPI::COMM_WORLD.Recv Recibe un mensaje MPI::Finalize Termina al MPI Estructura del Programa Maestro-Esclavo La estructura del programa se realizo para que el nodo maestro mande trabajos de manerasíncrona a los nodos esclavos. Cuando los nodos esclavos terminan la tarea asignada, avisan al nodomaestroparaqueseleasigneotratarea(estastareassonacordesalaetapa correspondiente del método de descomposición de dominio ejecutándose en un instantedado).enlamedidadeloposiblesetratademandarpaquetesdedatos acadanodoesclavoyqueestosregresentambiénpaquetesalnodomaestro,a manera de reducir las comunicaciones al mínimo y tratar de mantener siempre ocupados a los nodos esclavos para evitar los tiempos muertos, logrando con ello una granularidad gruesa, ideal para trabajar con clusters. La estructura básica del programa bajo el esquema maestro-esclavo codificadaenc++yusandompies: main(int argc, char*argv[]) { } MPI::Init(argc,argv); ME_id = MPI::COMM_WORLD.Get_rank(); MP_np = MPI::COMM_WORLD.Get_size(); if(me_id==0){ }else{ } // Operaciones del Maestro // Operaciones del esclavo con identificador ME_id MPI::Finalize(); En este único programa se deberá de codificar todas las tareas necesarias paraelnodomaestroycadaunodelosnodosesclavos,asícomolasformasde intercomunicación entre ellos usando como distintivo de los distintos procesos alavariableme_id.paramásdetallesdeestaformadeprogramaciónyotras funcionesdempiver[12]y[13]. Los factores limitantes para el esquema maestro-esclavo pueden ser de dos tipos, los inherentes al propio esquema maestro-esclavo y al método de descomposición de dominio: 8485 El esquema de paralelización maestro-esclavo presupone contar con un nodo maestro lo suficientemente poderoso para atender simultáneamente las tareas síncronas del método de descomposición dedominio,yaqueestedistribuyetareasacordealnúmerodesubdominios, estas si son balanceadas ocasionaran que todos los procesadores esclavos terminen al mismo tiempo y el nodo maestro tendrá que atender múltiples comunicaciones simultáneamente, degradando su rendimiento al aumentar el número de subdominios. Al ser síncrono el método de descomposición de dominio, si un nodoesclavoacabalatareaasignadayavisaalnodomaestro,este no podrá asignarle otra tarea hasta que todos los nodos esclavos concluyan la suya. Para los factores limitantes inherente al propio esquema maestro-esclavo, es posible implementar algunas operaciones del nodo maestro en paralelo, ya seausandoequiposmultiprocesadoroenmásdeunnododistintosalosnodos esclavos. Para la parte inherente al método de descomposición de dominio, la parte medularladaelbalanceodecargas.esdecirquecadanodoesclavotengauna cargadetrabajoigualalrestodelosnodos.estebalanceodecargaspuedeno ser homogéneo por dos razones: Al tener P procesadores en el equipo paralelo, la descomposición del dominio no sea la adecuada. Si se tiene una descomposición particular, esta se implemente en un número de procesadores inadecuado. Cualquiera de las dos razones generarán desbalanceo de la carga en los nodos esclavos, ocasionando una perdida de eficiencia en el procesamiento de un problema bajo una descomposición particular en una configuración del equipo paraleloespecifica,esporestoqueenalgunoscasosalaumentarelnúmerode procesadores que resuelvan la tarea no se aprecia una disminución del de tiempo de procesamiento. ElnúmerodeprocesadoresP queseusenpararesolverundominioytener buen balance de cargas puede ser conocido si aplicamos el siguiente procedimiento:sieldominiosedescomponeenn msubdominios(laparticióngruesa), entoncessegeneraráns=n msubdominios i,enestecaso,setieneunbuen balanceodecargassi(p 1) s.laparticiónfinaseobtienealdescomponera cadasubdominio i enp qsubdominios. Como ejemplo, supongamos que deseamos resolver el dominio usando 81 81nodos(nodos=n p+1ynodos=m q+1),demanerainmediatanossurgen las siguientes preguntas: cuales son las posibles descomposiciones validas? y en cuantos procesadores se pueden resolver cada descomposición?. Para este ejemplo, sin hacer la tabla exhaustiva obtenemos: 8586 Partición Subdominios Procesadores 1x2 y 80x40 2 2,3 1x4 y 80x20 5 2,5 1x5 y 80x16 6 2,6 2x1 y 40x80 2 2,3 2x2 y 40x40 4 2,3,5 2x4 y 40x20 8 2,3,5,9 2x5 y 40x ,3,6,11 2x8 y 40x ,3,5,9,17 4x1 y 20x80 4 2,3,5 4x2 y 20x40 8 2,3,5,9 4x4 y 20x ,3,5,9,17 4x5 y 20x ,3,5,6,11,21 5x1 y 16x80 5 2,6 5x2 y 16x ,3,6,11 5x4 y 16x ,3,5,6,11,21 5x5 y 16x ,6,26 De esta tabla es posible seleccionar (para este ejemplo en particular), las descomposiciones que se adecuen a las necesidades particulares del equipo con que se cuente. Sin embargo hay que tomar en cuenta siempre el número de nodos por subdominio de la partición fina, ya que un número de nodos muy grande puede que exceda la cantidad de memoria que tiene el nodo esclavo y un número pequeño estaría infrautilizando el poder computacional de los nodos esclavos. De las particiones seleccionadas se pueden hacer corridas de prueba para evaluar su rendimiento, hasta encontrar la que menor tiempo de ejecución consuma, maximizando así la eficiencia del equipo paralelo. Programación Paralela en Multihilos En una computadora, sea secuencial o paralela, para aprovechar las capacidades crecientes del procesador, el sistema operativo divide su tiempo de procesamiento entre los distintos procesos, de forma tal que para poder ejecutar a un proceso, el kernel les asigna a cada procesouna prioridad yconellounafraccióndel tiempototal de procesamiento,deformatalquesepuedaatenderatodosycadaunodelosprocesos de manera eficiente. En particular, en la programación en paralelo usando MPI, cada proceso(que eventualmente puede estar en distinto procesador) se lanza como una copia del programa con datos privados y un identificador del proceso único, de tal forma que cada proceso sólo puede compartir datos con otro proceso mediante paso de mensajes. Estaformadelanzarprocesosporcadatareaquesedeseehacerenparaleloes costosa, por llevar cada una de ellas todo una gama de subprocesos para poderle asignar recursos por parte del sistema operativo. Una forma más eficiente de hacerlo es que un proceso pueda generar bloques de subprocesos que puedan ser ejecutados como parte del proceso(como subtareas), así en el tiempo asignado 8687 se pueden atender a más de un subproceso de manera más eficiente, esto es conocido como programación multihilos. Los hilos realizarán las distintas tareas necesarias en un proceso. Para hacer que los procesos funcionen de esta manera, se utilizan distintas técnicas que le indican kernel cuales son las partes del proceso que pueden ejecutarse simultáneamente y el procesador asignará una fracción de tiempo exclusivo al hilo del tiempo total asignado al proceso. Los datos pertenecientes al proceso pasan a ser compartidos por los subprocesos lanzados en cada hilo y mediante una técnica de semáforos el kernel mantiene la integridad de estos. Esta técnica de programación puede ser muy eficiente si no se abusa de este recurso, permitiendo un nivel más de paralelización en cada procesador. Esta forma de paralelización no es exclusiva de equipos multiprocesadores o multicomputadoras, ya que pueden ser implementados a nivel de sistema operativo. 8788 8. Implementación Computacional Secuencial y Paralela de DDM Apartirdelosmodelosmatemáticos(capítulo2y3)ylosmodelosnuméricos(capítulos 4, 5, y 6), en este capítulo se describe el modelo computacional contenido en un programa de cómputo orientado a objetos en el lenguaje de programaciónc++ensuformasecuencialyensuformaparalelaenc++usando la interfaz de paso de mensajes(mpi) bajo el esquema maestro-esclavo (capítulo 7). Esto no sólo nos ayudará a demostrar que es factible la construcción del propio modelo computacional a partir del modelo matemático y numérico para la solución de problemas reales. Además, se mostrará los alcances y limitaciones en el consumo de los recursos computacionales, evaluando algunas de las variantes de los métodos numéricos con los que es posible implementar el modelo computacional y haremos el análisis de rendimiento sin llegar a ser exhaustivo esté. También exploraremos los alcances y limitaciones de cada uno de los métodos implementados(fem, DDM secuencial y paralelo) y como es posible optimizar los recursos computacionales con los que se cuente. Primeramente hay que destacar que el paradigma de programación orientada a objetos es un método de implementación de programas, organizados como colecciones cooperativas de objetos. Cada objeto representa una instancia de alguna clase y cada clase es miembro de una jerarquía de clases unidas mediante relaciones de herencia, contención, agregación o uso. Esto nos permite dividir en niveles la semántica de los sistemas complejos tratando así con las partes, que son más manejables que el todo, permitiendo su extensión y un mantenimiento más sencillo. Así, mediante la herencia, contención, agregación o usó nos permite generar clases especializadas que manejan eficientemente la complejidad del problema. La programación orientada a objetos organiza un programa entorno a sus datos (atributos) y a un conjunto de interfases bien definidas para manipular estos datos(métodos dentro de clases reusables) esto en oposición a los demás paradigmas de programación. El paradigma de programación orientada a objetos sin embargo sacrifica algo de eficiencia computacional por requerir mayor manejo de recursos computacionales al momento de la ejecución. Pero en contraste, permite mayor flexibilidad al adaptar los códigos a nuevas especificaciones. Adicionalmente, disminuye notoriamente el tiempo invertido en el mantenimiento y búsqueda de errores dentro del código. Esto tiene especial interés cuando se piensa en la cantidad de meses invertidos en la programación comparado con los segundos consumidos en la ejecución del mismo. Para empezar con la implementación computacional, primeramente definiremos el problema a trabajar. Este, pese a su sencillez, no pierde generalidad permitiendo que el modelo mostrado sea usado en muchos sistemas de la ingeniería y la ciencia. 8889 8.1. ElOperadordeLaplaceylaEcuacióndePoisson Consideramos como modelo matemático el problema de valor en la frontera (BVP) asociado con el operador de Laplace en dos dimensiones, el cual en general es usualmente referido como la ecuación de Poisson, con condiciones de frontera Dirichlet, definido en como: 2 u = f en (233) u = g en. Se toma está ecuación para facilitar la comprensión de las ideas básicas. Es un ejemplo muy sencillo, pero gobierna los modelos de muchos sistemas de la ingenieríaydelaciencia,entreelloselflujodeaguasubterráneaatravésdeun acuífero isotrópico, homogéneo bajo condiciones de equilibrio y es muy usada en múltiples ramas de la física. Por ejemplo, gobierna la ecuación de la conducción de calor en un sólido bajo condiciones de equilibrio. En particular consideramos el problema con definido en: =[ 1,1] [0,1] (234) donde cuya solución es f =2n 2 π 2 sin(nπx) sin(nπy) y g =0 (235) u(x, y) = sin(nπx) sin(nπy). (236) Paralaspruebasderendimientoenlascualesseevalúaeldesempeñodelos programasrealizadosseusan=10,peroesposiblehacerloconn Ngrande. Porejemploparan=4,lasoluciónesu(x,y)=sin(4πx) sin(4πy),cuyagráfica se muestra a continuación: Figura10:SoluciónalaecuacióndePoissonparan=4. 89 Mostrar más
RESUMEN DEL CURSO DE METODOS NUMERICOS impartido por Virginia Muto Foresi Departamento de Matemática Aplicada y Estadística e Investigación Operativa Facultad de Ciencia y Tecnología Universidad del País Más detalles Notas de CFD. Rev 0.2.2. Adrián Lozano Durán. adrian@torroja.dmt.upm.es
Notas de CFD Rev 0.2.2 Adrián Lozano Durán adrian@torroja.dmt.upm.es 24 de septiembre de 2013 Índice Índice 1 1 Computación Científica 4 1.1 El ordenador como herramienta para resolver problemas matemáticos Más detalles SIMULACIÓN. Texto Base de Teoría. Alfonso Urquía Moraleda
SIMULACIÓN Texto Base de Teoría Alfonso Urquía Moraleda Departamento de Informática y Automática Escuela Técnica Superior de Ingeniería Informática, UNED Juan del Rosal 16, 28040 Madrid, España E-mail: Más detalles Proyecto Fin de Carrera
Simulación numérica de la distribución de temperaturas y flujo de aire en el horno de una cocina de uso doméstico UNIVERSIDAD POLITÉCNICA SALESIANA SEDE CUENCA UNIDAD DE POSGRADOS MAESTRÍA EN MÉTODOS Más detalles Diseño y Construcción de Bases de Datos Distribuidas Heterogéneas sobre Oracle Y SQL Server
Universidad Carlos III de Madrid Escuela Politécnica Superior Ingeniería Técnica en Informática de Gestión Proyecto Fin de Carrera Diseño y Construcción de Bases de Datos Distribuidas Heterogéneas sobre Más detalles DATA WAREHOUSE PARA LA GESTIÓN DE LISTA DE ESPERA SANITARIA
UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO FACULTAD DE INGENIERÍA T E S I S I M P L E M E N T A C I Ó N D E U N S I S T E M A S E R V I C E D E S K B A S A D O E N I T I L. PARA OBTENER EL TÍTULO DE: INGENIERO Más detalles PONTIFICIA UNIVERSIDAD CATÓLICA DEL PERÚ FACULTAD DE CIENCIAS E INGENIERÍA
TESIS PUCP Esta obra ha sido publicada bajo la licencia Creative Commons Reconocimiento-No comercial-compartir bajo la misma licencia 2.5 Perú. Para ver una copia de dicha licencia, visite http://creativecommons.org/licenses/by-nc-sa/2.5/pe/ Más detalles Algoritmos evolutivos multiobjetivo para despacho de tareas en redes heterogéneas
Instituto de Computación Facultad de Ingeniería Universidad de la República Montevideo, Uruguay Proyecto de Grado Algoritmos evolutivos multiobjetivo para despacho de tareas en redes heterogéneas Darío Más detalles Universidad Nacional del Nordeste Facultad de Ciencias Exactas y Naturales y Agrimensura. Monografía de Adscripción: Data Warehouse
Universidad Nacional del Nordeste Facultad de Ciencias Exactas y Naturales y Agrimensura Monografía de Adscripción: Data Warehouse Rojas, Mariana Isabel LU: 38382 Prof. Director: Mgter. David Luis La Red Más detalles Introducción a los métodos de descomposición de dominio
Introducción a los métodos de descomposición de dominio Salim Meddahi Departamento de Matemáticas Universidad de Oviedo Francisco Javier Sayas Departamento de Matemática Aplicada Universidad de Zaragoza Más detalles Capítulo 1. Diseño de la LAN
1.0 Introducción del capítulo 1.0.1 Introducción del capítulo Para pequeñas y medianas empresas, la comunicación digital con datos, voz y video es esencial para la supervivencia de la empresa. En consecuencia, Más detalles DATA WAREHOUSE Y DATAMINING
CAPITULO IV DATA WAREHOUSE Y DATAMINING 4.1 Data Warehouse 4.2 DataMining Un concepto importante dentro de la integración de sistemas y, además, relacionado con los sistemas gerenciales es DataWarehouse. Más detalles Componentes software para gestión de dispositivos físicos. Estudio comparativo de soluciones basadas en.net y JavaBeans
UNIVERSIDAD POLITÉCNICA DE CARTAGENA E. T. S. Ingeniería de Telecomunicaciones Componentes software para gestión de dispositivos físicos. Estudio comparativo de soluciones basadas en.net y JavaBeans Andrés Más detalles UN PROCEDIMIENTO OPTIMAL PARA RESOLVER EL MEDIAN SHORTEST PATH PROBLEM
UNIVERSIDAD DEL BÍO-BÍO FACULTAD DE INGENIERÍA DEPARTAMENTO DE INGENIERÍA INDUSTRIAL PROFESOR GUÍA: SR. CARLOS OBREQUE NIÑEZ UN PROCEDIMIENTO OPTIMAL PARA RESOLVER EL MEDIAN SHORTEST PATH PROBLEM TRABAJO Más detalles INSTITUTO TECNOLOGICO SUPERIOR DE ZONGOLICA BASES DE DATOS DISTRIBUIDAS LSC. RAFAEL JUAN CORDOBA DE LA LLAVE
INSTITUTO TECNOLOGICO SUPERIOR DE ZONGOLICA BASES DE DATOS DISTRIBUIDAS LSC. RAFAEL JUAN CORDOBA DE LA LLAVE MAYO 2014 I N T R O D U C C I O N El mundo moderno se caracteriza por exceso de tecnología. Más detalles Francisco J. Ferrer Jesús Aguilar Joaquín Peña. Departamento de Lenguajes y Sistemas Informáticos Facultad de Informática Universidad de Sevilla
Francisco J. Ferrer Jesús Aguilar Joaquín Peña Departamento de Lenguajes y Sistemas Informáticos Facultad de Informática Universidad de Sevilla INFORME TÉCNICO: LSI-2000-08 Julio, 2000 1. Introducción Más detalles UNIVERSIDAD AUTÓNOMA DE BAJA CALIFORNIA SUR ÁREA DE CONOCIMIENTO DE CIENCIAS DEL MAR DEPARTAMENTO ACADÉMICO DE SISTEMASCOMPUTACIONALES TESIS
UNIVERSIDAD AUTÓNOMA DE BAJA CALIFORNIA SUR ÁREA DE CONOCIMIENTO DE CIENCIAS DEL MAR DEPARTAMENTO ACADÉMICO DE SISTEMASCOMPUTACIONALES TESIS CONSTRUCCIÓN DE UN DATAWAREHOUSE PARA DATOS DE PESCA RIBEREÑA Más detalles Arquitectura para un Sistema de Ficheros Distribuido Orientado al Ámbito de las Aplicaciones Grid y Cloud.
Arquitectura para un Sistema de Ficheros Distribuido Orientado al Ámbito de las Aplicaciones Grid y Cloud. Tesis de Máster (20 créditos). Máster Oficial: Computación Paralela y Distribuida. Autor: José Más detalles ESCUELA TÉCNICA SUPERIOR DE INGENIEROS INDUSTRIALES Y DE TELECOMUNICACIÓN MÁSTER UNIVERSITARIO EN TECNOLOGÍAS INFORMÁTICAS
ESCUELA TÉCNICA SUPERIOR DE INGENIEROS INDUSTRIALES Y DE TELECOMUNICACIÓN MÁSTER UNIVERSITARIO EN TECNOLOGÍAS INFORMÁTICAS Trabajo de Fin de Máster ESTUDIO DE MODELOS DE PREDICCIÓN DE CONSULTAS CONCURRENTES Más detalles Ingeniería de Tráfico en Redes MPLS
Ingeniería de Tráfico en Redes MPLS Proyecto Final de Carrera Integrantes: Adrián Delfino, Sebastián Rivero, Marcelo San Martín Tutor: Ing. Pablo Belzarena Instituto de Ingeniería Eléctrica Facultad de Más detalles UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS
UNIVERSIDAD DE CHILE FACULTAD DE CIENCIAS FÍSICAS Y MATEMÁTICAS DEPARTAMENTO DE CIENCIAS DE LA COMPUTACIÓN PARALELIZACIÓN DE UN ALGORITMO PARA LA DETECCIÓN DE CÚMULOS DE GALAXIAS MEMORIA PARA OPTAR AL Más detalles ANÁLISIS DE REQUISITOS
Índice 1 INTRODUCCIÓN...9 1.1 MOTIVACIÓN...9 1.1.1 Introducción...9 1.1.2 Razón y oportunidad...9 1.2 ESTUDIO DEL CASO...10 1.2.1 Historia y actualidad...10 1.2.2 En qué consiste el CRM?...12 1.2.3 Estudio Más detalles Trabajo de investigación I (Especialización Sistemas Distribuidos, Multimedia y Seguros) Autor: Eugenio Soler Galán Profesor: Benjamín Ramos Álvarez
http://www.eugenio-soler-galan.es/aleatoriedad Trabajo de investigación I Profesor: Benjamín Ramos Álvarez junio de 2008 Página 2 de 54 INDICE 1. MOTIVACIÓN... pág. 4 2. DÓNDE ENCONTRAMOS ALEATORIEDAD?... Más detalles INTRODUCCIÓN MARCO PROBLÉMICO JUSTIFICACIÓN 2. OBJETIVOS 2.1 OBJETIVO GENERAL 2.2 OBJETIVOS ESPECÍFICOS 3. MARCO REFERENCIAL 3.
INTRODUCCIÓN MARCO PROBLÉMICO JUSTIFICACIÓN 2. OBJETIVOS 2.1 OBJETIVO GENERAL 2.2 OBJETIVOS ESPECÍFICOS 3. MARCO REFERENCIAL 3.1 TÉCNICAS DE ALMACENAMIENTO DE DATOS 3.1.1 BASES DE DATOS 3.1.1.1 CARACTERÍSTICAS Más detalles ArchSync. Asistencia en la Sincronización de Documentación de Diseño con Implementación. por. Martín Blech y Juan Pablo Carlino
ArchSync Asistencia en la Sincronización de Documentación de Diseño con Implementación por Martín Blech y Juan Pablo Carlino Trabajo final de la carrera de grado de Ingeniero de Sistemas de la Universidad Más detalles Desarrollo de software dirigido por modelos para facilitar a usuarios inexpertos la aplicación de técnicas de minería de datos. Roberto Espinosa Oliva
Desarrollo de software dirigido por modelos para facilitar a usuarios inexpertos la aplicación de técnicas de minería de datos Roberto Espinosa Oliva Universidad de Alicante Instituto Universitario de Más detalles Análisis del Sistema de Información
Análisis del Sistema de Información ÍNDICE DESCRIPCIÓN Y OBJETIVOS... 2 ACTIVIDAD ASI 1: DEFINICIÓN DEL SISTEMA... 6 Tarea ASI 1.1: Determinación del Alcance del Sistema... 6 Tarea ASI 1.2: Identificación Más detalles Conjunto de cosas que ordenadamente relacionadas entre sí contribuyen a determinado Objeto
El concepto de sistema, término del cual tal vez se ha abusado en los últimos tiempos se aplica a los fenómenos más diversos y a veces sin demasiado rigor. Se trata de una noción difícil de precisar, posiblemente Más detalles 2016 © DocPlayer.es Política de privacidad | Condiciones del servicio | Feedback

References: resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución 
 resolución