Estudo comparativo de dois métodos econômicos aplicados a um esquema numérico de previsão do tempo

Comparative study of two economic methods applied to a numerical scheme of weather prediction

Estudo comparativo de dois métodos econômicos aplicados a um esquema numérico de previsão do tempo

A. Bourchtein

(CPMET/UFPEL)

Rua Anchieta 4715 bloco K, ap.304

96020-250 Pelotas-RS

[e-mail] Duas versões do esquema semi-implícito são apresentadas: o método implícito de direções alternadas (o método ADI) e a abordagem de Robert com o uso do algoritmo multigrade (o método MG). São realizados experimentos com o modelo regional baroclínico que usa 11 níveis para discretização vertical e passo 75 km para grades horizontais com 97x97 e 93x101 nós. A qualidade das previsões com 24 horas de antecedência não depende da versão do modelo. O método MG permite reduzir o tempo computacional gasto em torno de 14% em comparação com o método de Robert, que usa as relaxações sucessivas superiores. Para os mesmos passos temporais, o método multigrade de Robert é menos efetivo do que o método ADI. Mas o uso dos passos temporais máximos permitidos pela condição numérica da estabilidade de cada versão revela a vantagem da versão MG do esquema. A ausência da separação espacial permite usar os métodos multigrades no contexto de esquemas semi-implícitos semi-lagrangianos com passo temporal aumentado sem perda de precisão.

Palavras-chave: Modelo atmosférico baroclínico; Previsão numérica do tempo; Esquemas semi-implícitos; Métodos multigrades; Métodos implícitos de direções alternadas.

Comparative study of two economic methods applied to a numerical scheme of weather prediction - Two versions of the semi-implicit scheme using alternating directions implicit (ADI) method and Robert's approach with a multigrid (MG) method are presented. Experiments with a regional baroclinic model using 11 levels for the vertical discretization and 75 km step on horizontal grids with 97x97 and 93x101 points are carried out. The 24-hours forecasts quality does not depend on the model version. The MG method permits to reduce the computational time cost about 14% in comparison with Robert's method using successive over-relaxations. Under a choice of equal time steps Robert's MG method is less efficient than the ADI method. But usage of the maximum time steps allowed by the numerical stability condition of each version reveals an advantage of the MG version. Absence of the splitting permits to use the multigrid methods in frame of semi-implicit semi-lagrangian schemes with increased time step without loss of accuracy.

Key words: Atmospheric baroclinic model; Numerical weather prediction; Semi-implicit schemes; Multigrid methods; Alternating directions implicit methods.

INTRODUÇÃO

As duas principais exigências impostas a qualquer esquema numérico são a precisão na descrição dos fenômenos investigados e a eficiência computacional do esquema de integração. A precisão é fornecida pela resolução fina ("fine mesh model") e pela estabilidade do algoritmo numérico. A eficiência depende da complexidade do algoritmo usado e da presença de métodos da matemática numérica que permitam, de forma econômica, resolver os problemas numéricos surgidos em cada passo da integração no tempo.

O objetivo do modelo regional é a previsão das ondas sinóticas, isto é, das perturbações do estado básico da atmosfera cuja escala espacial é consideravelmente grande (comprimento característico horizontal de h₀ » 1000 km) e cujas mudanças são suficientemente lentas (período característico de t₀ » 24 h). Por isso as discretizações espacial e temporal serão concordantes se os passos t (pelo tempo) e h (pelo espaço) satisfizerem a seguinte relação aproximada:

Neste caso, os erros de truncamento espacial e temporal são comparáveis.

Geralmente é impossível escolher a relação indicada entre os passos t e h e ao mesmo tempo preservar a eficiência da passagem de um passo temporal ao seguinte. Isto é conseqüência da natureza das equações primitivas que sustentam as ondas gravitacionais rápidas, desprezíveis na previsão de escala sinótica, e da natureza dos problemas de diferenças finitas, os quais exigem determinada ligação entre t e h como conseqüência da condição da estabilidade linear do esquema. Formalmente esta ligação (normalmente é chamada condição de Courant-Friedrichs-Lewy - CFL) pode ser formulada pelo seguinte modo:

C £ 1,

onde é o número de Courant.

Praticamente para muitos esquemas aplicados às equações hidrostáticas o número de Courant pode ser escrito na forma

,

onde c_exp é a velocidade máxima das ondas usadas de modo explícito no esquema numérico.

As equações hidrostáticas linearizadas contém duas famílias de ondas mistas (gravitacionais-inerciais-advectivas) com velocidades até c_grav » 350 m/s e uma família de ondas advectivas com velocidades até c_adv» 80 m/s. Portanto a utilização dos esquemas semi-implícitos de Robert (eulerianos e semi-lagrangianos), que tratam de modo implícito os termos responsáveis pelas ondas mistas, é bem justificável, porque permite usar o passo temporal mais adequado ao passo espacial do que nos esquemas explícitos. Isso permite uma considerável diminuição na quantidade de passos temporais sem perda da qualidade dos prognósticos.

A parte implícita dos esquemas semi-implícitos em cada passo do tempo pode ser reduzida até uma equação elíptica tri-dimensional do tipo de Helmholtz. Usando os modos normais verticais do sistema discreto é possível reduzir o último problema ao conjunto das equações bi-dimensionais de Helmholtz. Certamente a eficiência dos esquemas depende essencialmente da existência de métodos econômicos da matemática numérica para resolução do problema linear bi-dimensional mencionado.

A utilização de métodos tradicionais como, por exemplo, SOR ("successive over-relaxation") gasta muito tempo de computador e diminui bastante a vantagem de esquemas semi-implícitos. Um dos métodos utilizados para desviar este problema é a combinação da técnica de separação com algoritmos semi-implícitos como nos artigos de Bates (1984), Tanguay & Robert (1986) e Bates & McDonald (1987). Este enfoque substitui o problema bi-dimensional pelo conjunto dos problemas uni-dimensionais resolvidos supereconomicamente através do método de varredura. O preço pago por essa eficiência é o erro adicional de truncamento que aumenta rapidamente com o aumento do passo t.

Outra abordagem efetiva contemporânea é a resolução direta dos problemas elípticos bi-dimensionais através dos métodos multigrades. Começando no ano de 1989, algumas variantes dos métodos multigrades foram utilizadas com sucesso em esquemas numéricos de previsão do tempo (Barros et al. 1989, Bates et al. 1990, Adams et al. 1992).

O objetivo deste trabalho é comparar a precisão e a eficiência dos métodos ADI e MG no contexto de um modelo semi-implícito baroclínico regional de previsão do tempo de curto prazo.

ESQUEMA SEMI-IMPLÍCITO DO MODELO BAROCLÍNICO

A aproximação semi-implícita de Robert et al. (1972) das equações hidrostáticas evolucionárias em coordenadas isobáricas numa projeção conforme do sistema esférico em sistema cartesiano é a seguinte:

Aqui , são as velocidades cartográficas; ; é um perfil médio da temperatura; ; as outras notações são as usuais. Para simplificar a forma das equações, é utilizado o fator de mapa igual a 1 e o relevo nulo.

O uso do esquema explícito "leap-frog"

permite reduzir o primeiro sistema ao sistema linear para as correções dq = q^t - q^e, q = u, v, f, T :

onde dq = q^e + q^-t -2q⁰ .

As equações diagnósticas da hidrostática e da continuidade são válidas em qualquer instante e por isso podem ser expressas na forma:

(2)

(3)

O resultado da discretização pela coordenada vertical p das equações (2) junto com a condição de contorno pode ser escrito na forma matricial:

onde é vetor-coluna de dimensão N dos valores da função q nos níveis prescritos no modelo (N é a quantidade dos níveis, N=11 no modelo atual) e B é matriz de ordem N .

Substituindo na equação do calor dT por dF através da equação da hidrostática (3), discretizando a equação obtida pela vertical e unindo com a condição na fronteira inferior (última equação em (1)), chegamos ao resultado:

. (5)

Excluindo de (5)

d

(6)

Unindo (6) com as primeiras duas equações de (1) discretizadas pela vertical, obtemos um sistema semi-discretizado de três equações para três incógnitas

(7)

Usando a expansão espectral da matriz C, C=S

S^-1SNCc_i²C

(8)

Aqui são coeficientes da expansão do vetor por autovetores da matriz C:

(9)

ou, como dizem freqüentemente, amplitudes dos modos normais verticais.

Depois da resolução dos sistemas (8) voltamos aos vetores , q = u, v, f usando as fórmulas (9). Dependendo do modo de resolução de cada sistema (8), por redução à equação bi-dimensional elíptica para ou por método ADI, obtemos uma versão do esquema semi-implícito de Robert ou esquema semi-implícito de direções alternadas.

MÉTODO ADI

Desde a decada de 70 o método ADI tem sido bastante usado em meteorologia. Suas diferentes versões são os elementos de vários esquemas apresentados nos artigos de Navon & Riphagen (1979), Bates (1984), Cohn et al. (1985) e outros.

A idéia do método é separar a parte implícita do sistema (8) por variáveis independentes x e y sem perda de aproximação e de estabilidade. Neste caso, a parte implícita em relação a cada variável pode ser reduzida à resolução do conjunto dos sistemas lineares com matrizes tri-diagonais, o que torna possível a aplicação do método econômico de varredura. As equações do método ADI são divididas em dois sistemas. O primeiro possui a parte implícita somente em relação à variável x (daqui por diante omitiremos o índice i e o símbolo _^ ) :

onde dq^p = q^p -q^e , q = u, v, f . Usando a 1ª equação excluímos du^p da 3ª e obtemos uma equação uni-dimensional em cada ponto fixo y para a variável df^p :

.

A aproximação espacial de 2ª ordem leva ao conjunto dos sistemas tri-diagonais para encontrar a função discreta incógnita dF^p.

O segundo sistema tem a forma:

onde dq^s = q^t - q^e , q = u, v, F . De modo análogo para o último sistema, em cada ponto fixo x, pode ser obtida a equação pela variável y:

cuja aproximação representa o conjunto dos sistemas tri-diagonais.

O esquema semi-implícito ADI é bastante efetivo, pois a quantidade de operações é diretamente proporcional à quantidade dos pontos da grade, e a análise linear da estabilidade numérica mostra que, utilizando um filtro de Asselin com o parâmetro

t

Aqui h é o passo horizontal, U é a velocidade máxima do vento.

MÉTODO DE ROBERT COM ALGORITMO MULTIGRADE

No método de Robert, usando a 1ª e 2ª equação do sistema (7), eliminamos as incógnitas

dudv

(12)

O esquema de Robert é estável linearmente se o passo t satisfaz a desigualdade:

(13)

A restrição (13) é mais suave do que a (11), mas o método de Robert exige, em cada passo do tempo, a solução do conjunto dos problemas bi-dimensionais de Helmholtz. As outras fórmulas do método são explícitas e resolvidas de modo simples. Por isso, a efetividade do método depende altamente da possibilidade de achar a solução do problema elíptico do tipo (12) de modo econômico.

Pelo que sabemos, atualmente não existe um método direto econômico de resolução da equação (12), isto é, um método cuja quantidade de operações é diretamente proporcional à quantidade de pontos da grade. Destes métodos os melhores são o método de transformação rápida de Fourier ("fast Fourier transformation") e a união deste com o método de varredura. Tais métodos exigem O(ML ln M) operações, onde MxL é a quantidade dos pontos da grade. Exceto os métodos MG, os melhores dos métodos iterativos (como ADI com passos efetivos, método de redução cíclico e outros) exigem o cálculo também com O(ML ln M) operações.

Os métodos MG foram propostos por Fedorenko (1962) e no mesmo artigo foi demonstrada sua economia teórica na resolução da equação de Poisson. O primeiro algoritmo econômico MG aplicado ao problema de Poisson foi apresentado no artigo de Brandt (1977).

Agora já existem diferentes programas dos métodos MG, onde a maioria dos quais é teoricamente econômica (Fulton et al. 1986). Cada programa contém várias versões e, para obter maior eficiência prática, é preciso levar em conta a especificidade do problema: tipo de operador elíptico, tipo de condições de contorno, dimensões da grade e a razão entre elas, tipo de aproximação inicial e outras. Portanto, há uma tarefa prática interessante de escolher uma variante do método MG mais efetiva para o problema elíptico dado.

Descrevemos brevemente a idéia dos métodos MG. Todos os métodos iterativos convergentes podem ser representados como um processo de suavização do erro inicial, isto é, a diminuição da diferença entre a solução exata do problema elíptico discreto e sua aproximação inicial.

Utilizando a análise de Fourier, pode-se mostrar que nesses processos a velocidade de convergência da solução é altamente dependente do comprimento de onda do erro inicial. Além disto, a velocidade de convergência mais lenta da "pior" onda (aquela com a menor taxa de convergência) é que determina a convergência média de todos os algoritmos.

A idéia principal dos métodos MG consiste em divisão de todo espectro ondular de erro inicial e posterior suavização de cada intervalo do espectro na grade respectiva, onde as ondas destacadas são mais curtas.

Um dos principais componentes do método MG é um algoritmo iterativo que suaviza as ondas curtas muito rapidamente. Se começarmos da grade mais fina (grade inicial com passo menor), após algumas iterações, as amplitudes dos erros iniciais das ondas curtas se tornam suficientemente pequenas. Para suavizar o próximo intervalo do espectro descemos à grade mais grossa (grade com maior passo entre os pontos), onde as ondas do intervalo escolhido são mais curtas. Por isso o mesmo algoritmo iterativo pode ser utilizado na grade grossa para a suavização das ondas mais longas. Este processo de aumento do passo das grades é repetido algumas vezes até que todas as partes do espectro sejam suavizadas, isto é, até descer à grade mais grossa possível. A volta para as grades mais finas é feita normalmente através de um dos métodos de interpolação.

Assim, cada parte do espectro do erro inicial é suavizada, na parte principal, na grade correspondente com a melhor velocidade de convergência. A quantidade de operações gastas em cada grade diminui rapidamente com o aumento do passo da grade. Portanto, economia computacional é obtida na utilização deste método.

RESULTADOS E DISCUSSÕES

Foi usado o pacote dos programas BOXMG, desenvolvido por Dendy (1982) e ajustado para os usuários por Bandy & Sweet (1992). O uso da discretização do tipo de Galerkin para a construção dos operadores em diferenças finitas para todas as grades, exceto a grade inicial, permite aplicar este pacote a problemas elípticos bi-dimensionais com qualquer número de nós.

Descrevemos algumas características do algoritmo usado:

1) na grade inicial foi escolhida a aproximação pelo esquema usual de 2ª ordem:

(14)

o algoritmo iterativo usado na grade inicial é relaxações bi-cores de Gauss-Seidel por pontos ("red-black point relaxations" - RB);

2) o pacote constrói automaticamente a seqüência das grades que assegura a efetividade do método MG independentemente do número de pontos;

3) em todas as grades com passo maior que h usam-se operadores aproximados com 9 pontos de padrão e relaxações quatro-cores de Gauss-Seidel;

4) aplica-se o esquema de correções ("Correction Scheme") no qual em cada grade mais grossa encontram-se as correções para a solução da grade mais fina;

5) o pacote oferece a possibilidade de escolher os parâmetros que determinam unicamente a versão do método MG. Os principais são: a escolha do método cíclico ("Cyclic MG") ou completo ("Full MG"); o tipo V ou W do ciclo MG ("V-cycle" ou "W-cycle"); o número dos ciclos; o número de relaxações numa grade; a aproximação inicial.

Através de experimentos, foi obtida uma variante efetiva que assegura a precisão exigida da solução com cálculo menor. Foi escolhido o método cíclico com uso do ciclo V. Foram aplicados dois ciclos MG para o 1º modo (com autovalor maior c₁» 300 m/s ) e um ciclo para todos os outros. A aproximação inicial nula foi usada para as primeiras horas da previsão e a correção do passo anterior foi aplicada como aproximação inicial para as últimas horas de previsão. Esta situação é ligada ao uso dos campos não inicializados como valores iniciais no modelo. Portanto ondas gravitacionais falsas com grande amplitude surgem no começo da previsão. Elas quase desaparecem depois de 15 horas devido às dissipações numérica e física usadas no modelo.

Mencionamos também, que o melhor resultado foi atingido usando um ciclo de relaxações em cada grade antes de passar para a grade mais fina ou mais grossa, exceto a grade fina inicial onde foram aplicadas duas relaxações.

As condições de fronteira para a equação (14) são as conseqüências da definição das funções na grade C de Arakawa que é usada no modelo. O deslocamento à direita da componente u e o deslocamento para baixo da componente v em relação ao geopotencial F leva às condições de Neumann (definição da derivada normal ¶F/¶n) nas fronteiras esquerda e superior. Nas fronteiras direita e inferior, usam-se as condições de Dirichlet (os próprios valores do geopotencial). A explicação deste fato encontra-se no Apêndice.

Na versão anterior do modelo para a resolução dos problemas de Helmholtz foi usado o método simples, usual e bastante rápido de relaxações sucessivas superiores ("successive overrelaxation" - SOR) que exige teoricamente O(ML²) operações. Para avaliar a vantagem do método MG em comparação com o método SOR foram escolhidos três pares das grades espaciais. O primeiro par é relacionado com o passo h=150 km e o número dos pontos MxL=49x49 (grade G₁ ) e MxL=47x51 (grade G₂ ). O segundo e terceiro par foram obtidos pela divisão sucessiva do passo espacial por 2, isto é, as grades G₃ e G₄ são relacionadas com o passo h=75 km e contém 97x97 e 93x101 nós respectivamente; as grades G₅e G₆ correspondem ao passo h=37,5 km e têm 193x193 e 185x201 nós.

O método tradicional para construção de grades grossas nos algoritmos MG é a duplicação do passo espacial. Neste caso existem restrições nas dimensões das grades: as vantagens dos algoritmos MG revelam-se somente quando ambas as dimensões M-1 e L-1são múltiplos de 2^m onde m é bastante grande. Por exemplo, usando a duplicação do passo inicial para a grade G₁ podem ser construídas mais quatro grades - a última, mais grossa, com resolução 4x4 pontos. Ao mesmo tempo a grade G₂ possui a única duplicação do passo - a grade mais grossa tem 24x26 nós. Portanto a aplicação do método tradicional para a grade G₁ é justificada, entretanto, para a grade G₂, não.

O pacote BOXMG usa uma discretização não tradicional do tipo de Galerkin o que permite escapar dessas restrições dimensionais. Comparações dos resultados de cálculos nas grades G₁ e G₂, G₃ e G₄ , G₅ e G₆ revelaram a independência prática da eficiência do algoritmo BOXMG em relação as dimensões da grade. Ao mesmo tempo, os experimentos nas grades G₁, G₂ e G₃ mostraram que o algoritmo BOXMG atinge eficiência não inferior do que os métodos multigrades tradicionais.

A comparação entre os três métodos (SOR, MG e ADI) foi realizada num modelo regional destinado para o Sul do Brasil. A região total de previsão ocupa aproximadamente uma área de 7000x7000 km². Ela está localizada, a princípio, entre 10 e 60 graus de latitudes meridionais e 10 e 100 graus de longuitudes oeste. Como dados iniciais, foram usados os campos de geopotencial e de vento de análise objetiva do NCEP USA. A região estimada de previsão fica dentro da região total com afastamento mínimo de 2000 km de fronteira da área total.

Uma avaliação da precisão do esquema foi realizada através de estimativas da qualidade dos prognósticos do geopotencial de 500 e 1000 mb. Com o objetivo de obter valores estatisticamente significativos, foi rodada uma série de 30 previsões com um prazo de antecedência de 24 horas para cada versão do modelo com passos h=150 km e h=75 km. Além disso foram rodadas duas previsões de cada versão com passo h=37,5 km. A Tab.1 mostra as correlações das tendências das previsões para passo h=75 km. As estimativas para o passo h=150 km são praticamente idênticas. Este resultado é conseqüência do uso da versão adiabática do modelo nos experimentos realizados.

Tabela 1 - Características dos esquemas com passo 75 km para previsões com prazo de antecedência de 24 horas. r500, r1000 são correlações das tendênciais para superfícies de 500 e de 1000 mb;t, K é passo temporal (seg) e número dos passos; CPU é tempo computacional gasto numa previsão (seg).Table 1 - 75 km schemes caracteristics for 24 hours integration. r500, r1000 are correlations of tendencies for 500 mb and 1000 mb surfaces; t, K is time step (sec) and number of steps; CPU is computer time cost for one forecast (sec).

As Fig.1 (para superfície de 500 mb) e Fig.2 (para superfície de 1000 mb) mostram os campos de análise objetiva, das diferenças entre análise e prognóstico e das diferenças entre versões distintas dos prognósticos para um dos dias da série. As diferenças máximas e médias quadráticas encontradas entre MG e SOR neste dia são 0,1 e 0,02 metros para ambas as superfícies (de 500 e de 1000 mb). As diferenças respectivas entre MG e ADI são 4,8 e 1,4 metros para a superfície de 500 mb e 6,1 e 1,5 metros para a de 1000 mb. A variação de dia a dia não superou 60% dos valores indicados. Assim as diferenças entre prognósticos podem ser consideradas desprezíveis.

Figura 1 - (a) Superfície de 500 mb de análise do NCEP para as 00:00 TMG do dia 30 de abril de 1996; (b) diferenças entre análise e prognóstico MG respectivo com 24 horas de antecedência; (c) diferenças entre prognósticos MG e SOR; (d) diferenças entre prognósticos MG e ADI. Os valores são dados em metros.

Figure 1 - (a) The 500 mb height NCEP analysis at 00:00 GMT 30 april 1996; (b) differences between the analysis and corresponding 24 hours MG forecast; (c) differences between MG and SOR forecasts; (d) differences between MG and ADI forecasts. The values are in meters.

Figura 2 - Mesmo que na

Fig.1

mas para superfície de 1000 mb.

Figura 2 - Same as Fig.1 but for the 1000 mb height.

Cada um dos três métodos (SOR, MG e ADI) foi aplicado na resolução de quatro problemas do tipo (12) com maiores autovalores c_i², pois os outros modos verticais são lentos (c_5-11<18m/s) e podem ser deixados sem correções. Nesta tarefa, nas grades com o passo h=150 km, o método BOXMG foi 1,8 vezes mais rápido do que o método SOR. Sob o uso das grades mais finas com passo h=75 km, o método BOXMG foi 2,1 vezes mais rápido do que o SOR. Tendo em vista que a resolução dos problemas (12) pelo método SOR ocupa aproximadamente 26% do tempo computacional do modelo baroclínico com resolução fina, é fácil prever que usando BOXMG em vez do SOR é possível economizar em torno de 14% do tempo computacional total. Realmente os experimentos com o modelo baroclínico mostraram que na grade fina a versão com SOR gasta 605 seg e a versão com BOXMG, 520 seg para a mesma qualidade de previsão. Se usarmos o passo espacial h=150 km, os valores respectivos são 63 seg e 60 seg. Finalmente, a mesma qualidade dos prognósticos foi atingida nas grades com passo h=37,5 km usando 77 min do tempo computacional no método MG e 94 min no método SOR (em torno de 18% de economia). Todos os experimentos foram realizados no computador DEC-3000.

A economia computacional conseguida com o uso de BOXMG não atingiu o resultado do modelo semi-implícito ADI, que gasta somente 56 seg nas grades com passo h=150 km e 455 seg nas grades finas com h=75 km. A comparação dos prognósticos nos métodos ADI e MG mostrou sua coincidência prática. Mas é preciso destacar que este resultado foi conseguido com o uso dos mesmos passos temporais t =720 seg para h=150 km e t =360 seg para h=75 km) em ambos os métodos. Entretanto, estes passos temporais são máximos para o modelo ADI, mas no modelo de Robert eles podem ser aumentados em vezes em correspondência aos critérios de estabilidade (11) e (13). Para a = 0,2 , usado no modelo, o coeficiente de aumento do passo é aproximadamente 1,22 . Isso significa que pode-se diminuir a quantidade de operações em 1,22 vezes. Então o tempo computacional gasto no modelo de Robert na grade fina é 425 seg e menor do que no modelo ADI. Os melhores parâmetros dos esquemas MG, SOR e ADI para o passo h=75 km são colocados na Tab.1. As tentativas de diminuir o parâmetro a no filtro de Asselin e usar o maior passo temporal no modelo ADI não foram bem sucedidas porque levaram a instabilidade numérica do esquema.

Levando em conta o resultado teórico de Tanguay & Robert (1986) confirmado praticamente em vários esquemas (Bates 1984, Cohn et al. 1985 e outros) sobre a ocorrência de queda rápida na precisão dos prognósticos no método ADI, quando o passo do tempo se eleva aproximadamente 30 minutos, o método MG pode ser considerado mais útil.

Finalmente salientamos que os resultados obtidos são válidos somente para o caso do modelo regional ou do modelo de mesoescala, em que o fator de mapa não varia consideravelmente. Caso contrário, surge a necessidade do uso de algum algoritmo MG mais sofisticado, que poderia levar até à perda de eficiência do cálculo. Um exemplo deste tipo foi apresentado no artigo de Bates et al. (1990), no qual a parte implícita do modelo barotrópico global foi reduzida até a equação de Helmholtz em toda a esfera com variação considerável dos coeficientes por causa da mudança crucial na relação entre os passos pela longitude e latitude. Isso levou à necessidade do uso de um procedimento de suavização com linhas-relaxações ("line-relaxations") e com direções diferentes nas regiões distintas. O algoritmo foi rejeitado posteriormente por causa da inflexibilidade dimensional e da eficiência não suficiente em favor de um algoritmo direto de transformação rápida de Fourier (Bates et al. 1993).

CONCLUSÃO

As equações baroclínicas adiabáticas foram integradas numa área limitada usando dois métodos econômicos (MG e ADI) no contexto do esquema semi-implícito. Ambos os algoritmos têm critério de estabilidade relacionado com a velocidade do vento e exigem MxL operações para cada passo temporal.

A precisão e eficiência dos dois métodos foi avaliada através de integrações com um prazo de antecedência de 24 horas realizadas na grade C com passos espaciais diferentes. Como dados iniciais foram usados os campos de geopotencial e vento de análise objetivo do NCEP USA. Uma variante mais efetiva do método MG foi escolhida. A comparação da última com o método SOR demonstrou a vantagem evidente do método MG. A comparação entre a versão mais efetiva dos métodos MG e o método ADI revelou certa vantagem da técnica de separação para o mesmo passo temporal. A escolha do seu próprio passo t permitido pela condição da estabilidade linear levou o método MG a ser um pouco mais efetivo.

Ambos os métodos em consideração podem ser incorporados de modo semelhante nos esquemas semi-implícitos semi-lagrangianos que possibilitam aumentar o passo t até 1 hora, pelo menos para a parte dinâmica do modelo. Pretendemos realizar essas versões do modelo futuramente e, particularmente, avaliar a influência do erro de separação na precisão dos cálculos.

APÊNDICE

Apresentamos a dedução das condições de fronteira na grade C para o problema elíptico (14) num caso simples do sistema:

(A1)

Discretizamos (A1) na grade C uni-dimensional usando o esquema semi-implícito:

(A2)

com condições nas extremidades:

(A3)

Substituindo as expressões para u_j^t da primeira equação (A2) e das condições (A3) na segunda equação (A2) obtemos o seguinte sistema para incógnitas F_j^t :

(A4)

onde

O sistema (A4) não inclui o valor F₁^t e para seu fechamento é suficiente usar a condição F_N t =H_N . Devido ao valor F₁^t ser desacoplado do sistema (A4) ele pode ser escolhido arbitrariamente. Por exemplo, adotando a condição F₁^t = F₂^t ou, o que é equivalente, a condição de Neumann em diferenças (F₁^t = F₂^t )/h=0 podemos reescrever o sistema (A4) na forma:

com condições complementares de Neumann à esquerda e de Dirichlet à direita:

Essa demonstração facilmente se estende ao caso bi-dimensional.

Submetido em: 08/01/98

Revisado pelo(s) autor(es) em: 25/05/98

Aceito em: 20/06/98

COMPARATIVE STUDY OF TWO ECONOMIC METHODS APPLIED TO A NUMERICAL SCHEME OF WEATHER PREDICTION

The use of the vertical normal modes of the discrete model allows to transform the 3-dimensional linear elliptic system for implicit part in a set of uncoupled 2-dimensional problems. Both the alternating directions and multigrid methods are applied to correct only the first four vertical modes of the model, because the other ones are slow and can be treated by the explicit mode, using the leap-frog scheme.

The application of the alternating directions implicit method to the linear terms allows transforming the implicit part of the scheme in a set of 1-dimensional elliptic problems. The last ones are solved by the economical sweeping algorithm (Thomas algorithm).

The approximation of the same linear terms by the Crank-Nicolson method leads to the set of 2-dimensional Helmholtz’s problems. Their solution in a grid of, approximately, 100x100 points by the classic successive overrelaxations algorithm spends much computational time. Multigrid methods offer the possibility to diminish the number of computation.

The time step allowed by the linear stability condition is a little bigger for the version of scheme with the use of multigrid methods than for the version with the use of the alternating directions implicit method. The prognostic fields are practically identical in both methods. So, the advantage of one of the versions depends essentially on the computational time cost for a prediction. Theoretically, the alternating directions method as the multigrid method are economic algorithms because their quantity of operations is directly proportional to number of points of the grid. But, in practice, the methods demand a different computational resource for fixed grids.

A kind of multigrid method was selected through numerical experiments as the most effective version for the elliptic problems described. The advantage of the multigrid method in comparison with the successive overrelaxation one is valued. The comparison between the multigrid and the alternating directions methods is performed for various space grid. For the fine mesh model in the context of the Eulerian semi-implicit scheme a quasi-equilibrium of both methods is demonstrated.

Both investigated algorithms can be included in a similar way in the semi-Lagrangian schemes which permit time steps up to 1 hour. In this case, as it was demonstrated in several works, the methods of separation in relation to time lose accuracy in consequence of the rapid increase of the truncation error when the time step is exceeded, approximately, 30 minutes. At the same time, multigrid methods assure an acceptable level of the truncation error for steps around 1 hour. Hence the last ones may present the preferred numerical technique for solving the elliptic problems related to the implicit part of the semi-Lagrangian semi-implicit models.

NOTA SOBRE AUTOR

NOTE ABOUT THE AUTHOR

Andrei Bourchtein

M.Sc. em Matemática Aplicada em 1984 pela Universidade de Vladivostok, Russia; Ph.D. em Matemática Aplicada em 1990 pelo Centro Hidrometeorológico da Russia, Moscou. Pesquisador e pesquisador sênior, em quatro anos, no Instituto Hidrometeorológico de Vladivostok. Professor adjunto, em dois anos, na Universidade Federal de Vladivostok. Desde 1995, professor adjunto nos cursos de pós-graduação em meteorologia na Universidade Federal de Pelotas. Áreas de pesquisa incluem meteorologia dinâmica, modelagem atmosférica e esquemas numéricos da hidrodinâmica.

• ADAMS, J., GARSIA, R., GROSS, B., HACK, J., HAIDVOGEL, D. & PIZZO, V. - 1992 - Applications of multigrid software in the atmospheric sciences. Mon. Wea. Rev., 120: 1447-1458.
• BANDY, V. & SWEET, R. - 1992 - A set of three drivers for BOXMG: a black box multigrid solver. Comm. Appl. Num. Methods, 8: 563-571.
• BARROS, S. R. M., DEE, D. P. & DICKSTEIN, F. - 1989 - A multigrid solver for semi-implicit global shallow water models. Atmos.-Ocean, 28: 24-47.
• BATES, J. R. - 1984 - An efficient semi-Lagrangian and alternating direction implicit method for integrating the shallow water equations. Mon. Wea. Rev., 112: 2033-2047.
• BATES, J. R. & MCDONALD, A. - 1987 - A semi-Lagrangian and alternating direction implicit method for integrating a multilevel primitive equation model. Short and Medium Range Numerical Weather Prediction. T.Matsuno, Universal Academy Press, Japan, 223-231.
• BATES, J. R., SEMAZZI, F. H. M., HIGGINS, R. W. & BARROS, S. R. M. - 1990 - Integration of the shallow water equations on the sphere using a vector semi-Lagrangian scheme with a multigrid solver. Mon. Wea. Rev., 118: 1615-1627.
• BATES, J. R., MOORTHI S. & HIGGINS R. W. - 1993 - A global multilevel atmospheric model using a vector semi-Lagrangian finite-difference scheme. Part I: Adiabatic formulation. Mon. Wea. Rew., 121: 244-263.
• BRANDT, A. - 1977 - Multi-level adaptive solutions to boundary value problems. Math. Comp., 31: 333-390.
• COHN, S. E., DEE, D., ISAACSON, E., MARCHESIN, D. & ZWAS, G. - 1985 - A fully implicit scheme for the barotropic primitive equations. Mon. Wea. Rev., 113: 436-448.
• DENDY, J. E., Jr. - 1982 - Black box multigrid. J. Comput. Phys., 48: 366-386.
• FEDORENKO, R. P. - 1962 - A relaxation method for solving elliptic diference equations. USSR Comput. Math. Math. Phys., 1: 1092-1096.
• FULTON, S. R., CIESIELSKI, P. E. & SCHUBERT, W. H. - 1986 - Multigrid methods for elliptic problems: a review. Mon. Wea. Rev., 114: 943-959.
• NAVON, I. M. & RIPHAGEN, H. A. - 1979 - An implicit compact fourth-order algorithm for solving the shallow-water equations in conservation-law form. Mon. Wea. Rev., 107: 1107-1127
• ROBERT, A., HENDERSON, J. & TURNBULL, C. - 1972 - An implicit time integration scheme for baroclinic models of the atmosphere. Mon. Wea. Rev., 100: 329-335.
• TANGUAY M. & ROBERT A. - 1986 - Elimination of the Helmholtz equation associated with the semi-implicit scheme in a grid point model of the shallow water equations. Mon. Wea. Rev., 114: 2154-2162.

Datas de Publicação

• Publicação nesta coleção
  23 Dez 1999
• Data do Fascículo
  Mar 1998

Histórico

• Aceito
  20 Jun 1998
• Revisado
  25 Maio 1998
• Recebido
  08 Jan 1998