Métodos de resolução de sistemas não-lineares de difusão por advecção além de Newton-Raphson?

9

Estou trabalhando em um projeto no qual tenho dois domínios acoplados adv-diff por meio de seus respectivos termos de origem (um domínio adiciona massa, o outro subtrai massa). Por uma questão de brevidade, estou modelando-os em estado estacionário. As equações são sua equação de transporte de difusão de advecção padrão com um termo de origem semelhante a este:

c1t=0=F1+Q1(c1,c2)c2t=0=F2+Q2(c1,c2)

Onde é fluxo difusivo e advective para espécies , e é o termo fonte de espécies .FiiQii

Consegui escrever um solucionador para o meu problema usando o método Newton-Raphson e acoplei completamente os dois domínios usando uma matriz de massa de bloco, ou seja:

Fcoupled=[A100A2][c1,ic2,i]xi[b1(c1,i,c2,i)b2(c1,i,c2,i)]

O termo é usado para determinar a matriz jacobiana e atualizar e :Fcoupledc1c2

J(xi)[xi+1xi]=Fcoupled

ou

xi+1=xi(J(xi))1Fcoupled

Para acelerar as coisas, eu não calculo o Jacobiano a cada iteração - agora estou jogando com cada cinco iterações, o que parece funcionar bem o suficiente e manter a solução estável.

O problema é: eu vou mudar para um sistema maior, onde os dois domínios estão em 2D / 2.5D, e o cálculo da matriz jacobiana esgotará rapidamente meus recursos de computador disponíveis. Estou construindo este modelo para ser usado em uma configuração de otimização mais tarde, então também não posso estar ao volante a cada iteração que ajusta o fator de amortecimento, etc.

Estou certo em procurar outro algoritmo mais robusto para o meu problema ou isso é o melhor possível? Analisei um pouco a quase-linearização, mas não tenho certeza de como isso é aplicável ao meu sistema.

Existem outros algoritmos lisos que eu posso ter esquecido que podem resolver um sistema de equações não lineares sem recorrer a recalcular o jacobiano como ofensivo?

cbcoutinho
fonte
2
Você já considerou solucionadores iterativos, como métodos multigrid algébricos AMG. Pode ser necessário criar bons pré-condicionadores baseados na física.
precisa saber é o seguinte
11
Você pode obter acesso a um cluster de computação onde pode distribuir a formação e a solução jacobiana usando um pacote de álgebra linear paralela?
Bill Barth
Não, não considerei a AMG, pensei que eram apenas para sistemas simétricos e não podiam ser usados ​​em problemas dominados por convecção. Vou procurar novamente na literatura a AMG.
Cbcoutinho 27/04
Cálculos paralelos são difíceis porque este projeto está sendo desenvolvido como um aplicativo independente para colegas que não têm acesso a esse tipo de recursos. Eu considerei desenvolvimento MPI no projeto para o meu próprio bem, mas que iria aumentar a barreira de entrada para os outros, que foi o ponto inteiro em primeiro lugar ..
cbcoutinho
3
Por que a computação do jacobiano é tão problemática? Se você estiver usando diferenças / volumes / elementos finitos, ele deve ter uma parte esparsa que é sempre a mesma e uma parte diagonal que muda, mas é trivial de calcular.
David Ketcheson

Respostas:

4

Estou assumindo que a limitação em 2D e 3D é armazenar o jacobiano.

Uma opção é reter as derivadas de tempo e usar um "pseudo" passo explícito para iterar para o estado estacionário. Normalmente, o número de CFL necessário para sistemas difusivos e reativos pode ficar proibitivamente pequeno. Você pode tentar multigrid não linear (também chamado multigrid de armazenamento aproximado total) e tempo local para acelerar a convergência.

DF(un)δun=F(un)
DF
DF(un)δuF(un+ϵδuδu)F(un)ϵ.
AAxx

Agora, com um valor adequado de (geralmente cerca de para flutuadores de precisão dupla), é possível executar um loop de Newton sem nunca computar ou armazenar um jacobiano. Sei que essa técnica é usada para resolver alguns casos não triviais na dinâmica de fluidos computacional. Observe, no entanto, que o número de avaliações da função será maior do que em uma técnica de armazenamento em matriz, em vez de exigir um produto vetor de matriz.ϵ107F

Outra coisa a ser observada é que, se o seu sistema for necessário para um pré-condicionador poderoso (por exemplo, Jacobi ou Jacobi em bloco não será suficiente), convém tentar usar o método mencionado acima como mais suave em um esquema multigrid. Se você quiser experimentar um pré-condicionador Jacobi de ponto ou bloco, poderá calcular e armazenar apenas os elementos diagonais ou blocos diagonais do Jacobiano, o que não é muito. Eu também mencionaria que um pré-condicionador Gauss-Seidel ou SSOR pode ser possível implementar sem armazenar explicitamente um jacobiano. Este artigo descreve a implementação de um GMRES sem matriz pré-condicionado com Gauss-Seidel simétrico sem matriz no contexto da dinâmica de fluidos computacional.

Aditya Kashi
fonte
1

Da minha experiência com as equações de Navier-Stokes, pode-se fazer muito bem sem esquemas totalmente implícitos.

Se você deseja apenas um esquema numérico rápido para a solução da evolução temporal, dê uma olhada nos esquemas IMEX (implícito-explícito); veja, por exemplo, este artigo de Ascher, Ruuth, Spiteri Métodos Runge-Kutta implícito-explícito para equações diferenciais parciais dependentes do tempo .

Você pode até tentar usar apenas um esquema explícito de integração de tempo de alta ordem com controle de tamanho de etapas (como o Matlab ODE45). No entanto, você pode ter problemas devido à rigidez do seu sistema, proveniente da parte difusiva. Felizmente, a parte difusiva é linear para que possa ser tratada implicitamente (que é a idéia dos esquemas IMEX).

Jan
fonte
0

Inicialmente, considerei adicionar apenas uma observação, mas o espaço não era suficiente, por isso acrescento uma breve descrição de minhas experiências com esse tópico.

Primeiro, olhando para sua notação de , não vejo a forma acoplada, suponho que e sejam ambos dependentes de e . Além disso, se e são representações matriciais de aproximações de e então eles não devem depender apenas de , mas também de valores vizinhos, mas isso pode ser apenas um entendimento errado da sua notação .Fcoupledb1b2c1,ic2,iA1A2F1F2ci

Como comentário geral, gostaria de acrescentar que o uso de Jacobiano analítico parece ser a única maneira de obter uma convergência quadrática de solucionador iterativo não linear (por exemplo, o solucionador de Newton-Raphson no seu caso). Você observou isso no seu caso? É muito importante, porque, caso contrário, pode haver algum mal-entendido em suas aproximações (linearização).

Em todas as aplicações em que estive envolvido (algumas delas incluíam cálculos em larga escala), nunca tivemos problemas com o consumo de tempo ao montar o jacobiano, a questão mais demorada era sempre aplicar um solucionador linear. O jacobiano analítico (se disponível) sempre esteve em aplicações nas quais eu estava trabalhando na escolha preferida devido à convergência quadrática. Em alguns casos, esse solucionador não linear produz uma matriz que causa problemas na convergência do solucionador linear iterativo; portanto, tentamos usar uma linearização mais simples que a Jacobiana analítica para ajudar o solucionador linear. Essa troca entre o comportamento do solucionador algébrico não linear e linear, dependendo da linearização do sistema algébrico não linear, sempre foi complicada e não pude dar uma recomendação geral.

Mas você está certo de que a desvantagem (ou propriedade) do jacobiano analítico para o sistema de PDEs é que ele produz um sistema algébrico acoplado; portanto, se você dissocia esse sistema de alguma forma, por exemplo, resolvendo separadamente a aproximação de cada PDE por, digamos, divisão iterativa método, novamente você perde a convergência quadrática do solucionador global. Mas, pelo menos, se você resolver separadamente cada PDE discreto (dissociado), poderá acelerar novamente a solução para esse problema específico usando o método Newton-Raphson.

Peter Frolkovič
fonte
Olá, Peter, você está certo sobre o acoplamento. Editei a equação principal para mostrar e e são funções de e . As matrizes e , neste caso, são as matrizes de rigidez dos dois sistemas, desenvolvidas usando o método dos elementos finitos. Essas são apenas funções das coordenadas dos nós, e não das variáveis ​​de estado. e são vetores, portanto, são funções de um vetor de variáveis ​​de estado, não apenas de uma variável. Calculo numericamente um jacobiano usando diferenças finitas. Até agora, não investiguei um jacobiano analítico. b1b2c1c2A1A2F1F2
cbcoutinho