Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

| | | |

Matemática Numérica Paralela

3 Computação paralela e distribuída (MPI)3.4 Reduções

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

3.5 Aplicação: Equação do calor

Em revisão

Nesta seção, vamos construir um código MPI para a resolução da equação do calor pelo método das diferenças finitas. Como um caso teste, vamos considerar a equação do calor

u_{t}-u_{xx}=0,\quad 0<x<1,t>0,

(3.15)

com condições de contorno

u(0,t)=u(1,t)=0,t>0,

(3.16)

e condições iniciais

u(x,0)=\operatorname{sen}(\pi x),0\leq x\leq 1,

(3.17)

onde, $u=u(t,x)$ .

Seguindo o método de Rothe³³endnote: ³Erich Hans Rothe, 1895 - 1988, matemático alemão. Fonte: Wikipédia.⁴⁴endnote: ⁴Também chamado de método das linhas., vamos começar pela discretização no tempo. Para tanto, vamos usar a notação $t_{m}=m\cdot h_{t}$ , $m=0,1,2,\dotsc,M$ , onde $h_{t}$ é o passo no tempo e $M$ o número total de iterações no tempo. Usando o método de Euler⁵⁵endnote: ⁵Leonhard Paul Euler, 1707 - 1783, matemático suíço. Fonte: Wikipédia., a equação (3.15) é aproximada por

u(t_{m+1},x)\approx u(t_{m},x)+h_{t}u_{xx}(t_{m},x),\quad 0<x<1,

(3.18)

onde $u(t_{0},x)=u(0,x)$ , dada pela condição inicial (3.17).

Agora, vamos denotar $x_{i}=ih_{x}$ , $i=0,1,2,\dotsc,I$ , onde $h_{x}=1/I$ é o tamanho da malha (passo espacial). Então, usando a aproximação por diferenças finitas centrais de ordem 2, obtemos

u(t_{m+1},x_{i})\approx u(t_{m},x_{i})+h_{t}\left[\frac{u(t_{m},x_{i-1})-2u(t_% {m},x_{i})+u(t_{m},x_{i+1})}{h_{x}^{2}}\right],

(3.19)

para $i=1,2,\dotsc,I-1$ , observando que, das condições de contorno (3.16), temos

u(t_{m},x_{0})=u(t_{m},x_{I})=0.

(3.20)

Em síntese, denotando $u^{m}_{i}\approx u(t_{m},x_{i})$ , temos que uma aproximação da solução do problema de calor acima pode ser calculada pela seguinte iteração

u^{m+1}_{i}=u^{m}_{i}+\frac{h_{t}}{h_{x}^{2}}u^{m}_{i-1}-2\frac{h_{t}}{h_{x}^{% 2}}u^{m}_{i}+\frac{h_{t}}{h_{x}^{2}}u^{m}_{i+1},

(3.21)

com $u^{0}_{i}=\operatorname{sen}(\pi x_{i})$ e $u^{m}_{0}=u^{m}_{I}=0$ .

Para construirmos um código MPI, vamos fazer a distribuição do processamento pela malha espacial entre os $n_{p}$ processos disponíveis. Mais precisamente, cada processo $p=0,1,2,\dotsc,n_{p}-1$ , computará a solução $u^{m}_{i}$ nos ponto de malha $i=i_{p},i_{p}+1,\dotsc,f_{p}$ , onde

	$\displaystyle i_{p}=p\left\lfloor\frac{I}{n_{p}}\right\rfloor,\quad p=0,1,% \dotsc,n_{p}-1,$		(3.22)
	$\displaystyle f_{p}=(p+1)\left\lfloor\frac{I}{n_{p}}\right\rfloor,\quad p=0,1,% \dotsc,n_{p}-2,$		(3.23)

com $f_{n_{p}-1}=I$ . Ainda, por simplicidade e economia de memória computacional, vamos remover o superíndice, denotando por $u$ a solução na iteração corrente e por $u^{0}$ a solução na iteração anterior.

Neste caso, a cada iteração $m=0,1,2,\ldots,M$ , o processo $p=0$ , irá computar

u_{j}=u^{0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j-1}-2\frac{h_{t}}{h_{x}^{2}}u^{% 0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j+1},

(3.24)

para $j=1,\dotsc,f_{0}$ , com as condições de contorno

	$\displaystyle u^{0}_{0}=0$		(3.25)
	$\displaystyle u^{0}_{f_{0}+1}=u^{0}_{i_{1}+1}.$		(3.26)

Ou seja, este passo requer a comunicação entre o processo 0 e o processo 1.

Os processos $p=1,2,\dotsc,n_{p}-2$ irão computar

u_{j}=u^{0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j-1}-2\frac{h_{t}}{h_{x}^{2}}u^{% 0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j+1},

(3.27)

para $j=i_{p}+1,\dotsc,f_{p}$ , com as condições de contorno

	$\displaystyle u^{0}_{i_{p}}=u^{0}_{f_{p-1}}$		(3.28)
	$\displaystyle u^{0}_{f_{p}+1}=u^{0}_{i_{p+1}+1}.$		(3.29)

Logo, este passo requer a comunicação entre os processos $p-1$ , $p$ e $p+1$ .

Ainda, o processo $p=n_{p}-1$ deve computar

u_{j}=u^{0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j-1}-2\frac{h_{t}}{h_{x}^{2}}u^{% 0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j+1},

(3.30)

para $j=i_{n_{p}-1}+1,\dotsc,f_{n_{p}-1}-1$ , com as condições de contorno

	$\displaystyle u^{0}_{i_{n_{p}-1}}=u^{0}_{f_{n_{p}-2}}$		(3.31)
	$\displaystyle u^{0}_{f_{n_{p}-1}}=0.$		(3.32)

O que requer a comunicação entre os processos $n_{p}-2$ e o $n_{p}-1$ .

O código abaixo calor.cc é uma implementação MPI deste algoritmo. Cada processo aloca e computa a solução em apenas um pedaço da malha, conforme descrito acima. As comunicações entre os processos ocorrem apenas nas linhas 60-80. Verifique!

Código: calor.cc

⬇

1#include <stdio.h>

2#include <math.h>

4// API MPI

5#include <mpi.h>

7int main (int argc, char** argv) {

8 // Inicializa o MPI

9 MPI_Init(NULL, NULL);

11 // numero total de processos

12 int world_size;

13 MPI_Comm_size(MPI_COMM_WORLD, &world_size);

15 // ID (rank) do processo

16 int world_rank;

17 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

19 // parametros

20 size_t M = 1000;

21 size_t I = 10;

23 // tamanho dos passos discretos

24 double ht = 1e-3;

25 double hx = 1.0/I;

26 double cfl = ht/(hx*hx);

28 // malha espacial local

29 size_t ip = world_rank * int (I/world_size);

30 size_t fp = (world_rank+1) * int (I/world_size);

31 if (world_rank == world_size-1)

32 fp = I;

33 size_t my_I = fp - ip;

35 double x[my_I+1];

36 for (size_t j=0; j<=my_I; j++)

37 x[j] = (ip+j) * hx;

39 // solucao local

40 double u0[my_I+1], u[my_I+1];

42 // condicao inicial

43 for (size_t j=0; j<=my_I; j++) {

44 u0[j] = sin (M_PI * x[j]);

45 }

47 // condicoes de contorno de Dirichlet

48 if (world_rank == 0)

49 u[0] = 0.0;

50 if (world_rank == world_size-1)

51 u[my_I] = 0.0;

53 // auxiliares

54 double u00, u0I;

57 // iteracoes no tempo

58 for (size_t m=0; m<M; m++) {

60 if (world_rank == 0) {

61 MPI_Send (&u0[my_I], 1, MPI_DOUBLE,

62 1, 0, MPI_COMM_WORLD);

63 MPI_Recv (&u0I, 1, MPI_DOUBLE,

64 1, 0, MPI_COMM_WORLD,

65 MPI_STATUS_IGNORE);

66 } else if (world_rank < world_size-1) {

67 MPI_Recv (&u00, 1, MPI_DOUBLE,

68 world_rank-1, 0, MPI_COMM_WORLD,

69 MPI_STATUS_IGNORE);

70 MPI_Send (&u0[my_I], 1, MPI_DOUBLE,

71 world_rank+1, 0, MPI_COMM_WORLD);

73 MPI_Recv (&u0I, 1, MPI_DOUBLE,

74 world_rank+1, 0, MPI_COMM_WORLD,

75 MPI_STATUS_IGNORE);

76 MPI_Send (&u0[1], 1, MPI_DOUBLE,

77 world_rank-1, 0, MPI_COMM_WORLD);

78 } else {

79 MPI_Recv (&u00, 1, MPI_DOUBLE,

80 world_size-2, 0, MPI_COMM_WORLD,

81 MPI_STATUS_IGNORE);

82 MPI_Send (&u0[1], 1, MPI_DOUBLE,

83 world_size-2, 0, MPI_COMM_WORLD);

84 }

86 // atualizacao

87 u[1] = u0[1]

88 + cfl * u00

89 - 2*cfl * u0[1]

90 + cfl * u0[2];

91 for (size_t j=2; j<my_I; j++)

92 u[j] = u0[j]

93 + cfl * u0[j-1]

94 - 2*cfl * u0[j]

95 + cfl * u0[j+1];

96 if (world_rank < world_size-1)

97 u[my_I] = u0[my_I]

98 + cfl * u0[my_I-1]

99 - 2*cfl * u0[my_I]

100 + cfl * u0I;

101

102 // prepara nova iteracao

103 for (size_t j=0; j<=my_I; j++)

104 u0[j] = u[j];

105

106 }

107

108

109 // Finaliza o MPI

110 MPI_Finalize ();

111

112 return 0;

113}

No código acima, as comunicações entre os processos foram implementadas de forma cíclica. A primeira a ocorrer é o envio de $u^{0}_{f_{p}}$ do processo zero (linhas 61-62) e o recebimento de $u^{0}_{i_{p}}$ pelo processo 1 (linhas 67-69). Enquanto essa comunicação não for completada, todos os processos estarão aguardando, sincronizados nas linhas 67-69 e 79-81. Ocorrida a comunicação entre os processos 0 e 1, ocorrerá a comunicação entre o 1 e o 2, entre o 2 e o 3 e assim por diante, de forma cíclica (linhas 67-71 e 79-81).

As últimas comunicações também ocorrem de forma cíclica, começando pelo envio de $u^{0}_{i_{p}+1}$ pelo processo $n_{p}-1$ (linhas 81-83) e o recebimento de $u^{0}_{f_{p}+1}$ pelo processo $n_{p}-2$ (linhas 73-75). Em sequência, ocorrem as comunicações entre o processo $n_{p}-2$ e o processo $n_{p}-3$ , entre o $n_{p}-3$ e o $n_{p}-4$ , até o término das comunicações entre o processo 1 (linhas 76-77) e o processo 0 (linhas 63-65).

3.5.1 Exercícios

Em revisão

E. 3.5.1.

O problema (3.15)-(3.17) tem solução analítica

u(t,x)=e^{-\pi^{2}t}\operatorname{sen}(\pi x).

(3.33)

Modifique o código calor.cc acima, de forma que a acada iteração $m$ no tempo, seja impresso na tela o valor de $u^{m}(0.5)$ computado e o valor da solução analítica $u(t_{m},0.5)$ . Faça refinamentos nas malhas temporal e espacial, até conseguir uma aproximação com três dígitos significativos corretos.

E. 3.5.2.

Modifique o código calor.cc acima, de forma que a cada iteração $m$ no tempo, seja impresso na tela o valor do funcional

q(t)=\int_{0}^{1}u(t,x)\,dx.

(3.34)

Valide os resultados por comparação com a solução analítica (veja o Exercício (3.5.1)).

E. 3.5.3.

Modifique o código calor.cc acima, de forma que as comunicações iniciem-se entre os processos $n_{p}-1$ e $n_{p}-2$ . Ou seja, que a primeira comunicação seja o envio de $u^{0}_{i_{p}+1}$ pelo processo $n_{p}-1$ e o recebimento de $u^{0}_{f_{p}+1}$ pelo processo $n_{p}-2$ . E que, de forma cíclica, as demais comunicações ocorram. Poderia haver alguma vantagem em fazer esta modificação?

E. 3.5.4.

Modifique o código calor.cc acima, de forma que as comunicações sejam feitas de forma assíncrona, i.e. usando as rotinas MPI_Isend e MPI_Irecv. Há vantagens em se fazer esta modificação?

E. 3.5.5.

Faça um código MPI para computar uma aproximação da solução de

u_{t}-u_{xx}=0,\quad 0<x<1,t>0,

(3.35)

com condições de contorno periódicas

u(0,t)=u(1,t),t>0,

(3.36)

e condições iniciais

u(x,0)=\operatorname{sen}(\pi x),0\leq x\leq 1,

(3.37)

onde, $u=u(t,x)$ .

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Matemática Numérica Paralela

3 Computação paralela e distribuída (MPI)3.4 Reduções

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

3.5 Aplicação: Equação do calor

Em revisão

Nesta seção, vamos construir um código MPI para a resolução da equação do calor pelo método das diferenças finitas. Como um caso teste, vamos considerar a equação do calor

u_{t}-u_{xx}=0,\quad 0<x<1,t>0,

(3.15)

com condições de contorno

u(0,t)=u(1,t)=0,t>0,

(3.16)

e condições iniciais

u(x,0)=\operatorname{sen}(\pi x),0\leq x\leq 1,

(3.17)

onde, $u=u(t,x)$ .

u(t_{m+1},x)\approx u(t_{m},x)+h_{t}u_{xx}(t_{m},x),\quad 0<x<1,

(3.18)

onde $u(t_{0},x)=u(0,x)$ , dada pela condição inicial (3.17).

Agora, vamos denotar $x_{i}=ih_{x}$ , $i=0,1,2,\dotsc,I$ , onde $h_{x}=1/I$ é o tamanho da malha (passo espacial). Então, usando a aproximação por diferenças finitas centrais de ordem 2, obtemos

u(t_{m+1},x_{i})\approx u(t_{m},x_{i})+h_{t}\left[\frac{u(t_{m},x_{i-1})-2u(t_% {m},x_{i})+u(t_{m},x_{i+1})}{h_{x}^{2}}\right],

(3.19)

para $i=1,2,\dotsc,I-1$ , observando que, das condições de contorno (3.16), temos

u(t_{m},x_{0})=u(t_{m},x_{I})=0.

(3.20)

Em síntese, denotando $u^{m}_{i}\approx u(t_{m},x_{i})$ , temos que uma aproximação da solução do problema de calor acima pode ser calculada pela seguinte iteração

u^{m+1}_{i}=u^{m}_{i}+\frac{h_{t}}{h_{x}^{2}}u^{m}_{i-1}-2\frac{h_{t}}{h_{x}^{% 2}}u^{m}_{i}+\frac{h_{t}}{h_{x}^{2}}u^{m}_{i+1},

(3.21)

com $u^{0}_{i}=\operatorname{sen}(\pi x_{i})$ e $u^{m}_{0}=u^{m}_{I}=0$ .

	$\displaystyle i_{p}=p\left\lfloor\frac{I}{n_{p}}\right\rfloor,\quad p=0,1,% \dotsc,n_{p}-1,$		(3.22)
	$\displaystyle f_{p}=(p+1)\left\lfloor\frac{I}{n_{p}}\right\rfloor,\quad p=0,1,% \dotsc,n_{p}-2,$		(3.23)

Neste caso, a cada iteração $m=0,1,2,\ldots,M$ , o processo $p=0$ , irá computar

u_{j}=u^{0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j-1}-2\frac{h_{t}}{h_{x}^{2}}u^{% 0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j+1},

(3.24)

para $j=1,\dotsc,f_{0}$ , com as condições de contorno

	$\displaystyle u^{0}_{0}=0$		(3.25)
	$\displaystyle u^{0}_{f_{0}+1}=u^{0}_{i_{1}+1}.$		(3.26)

Ou seja, este passo requer a comunicação entre o processo 0 e o processo 1.

Os processos $p=1,2,\dotsc,n_{p}-2$ irão computar

u_{j}=u^{0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j-1}-2\frac{h_{t}}{h_{x}^{2}}u^{% 0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j+1},

(3.27)

para $j=i_{p}+1,\dotsc,f_{p}$ , com as condições de contorno

	$\displaystyle u^{0}_{i_{p}}=u^{0}_{f_{p-1}}$		(3.28)
	$\displaystyle u^{0}_{f_{p}+1}=u^{0}_{i_{p+1}+1}.$		(3.29)

Logo, este passo requer a comunicação entre os processos $p-1$ , $p$ e $p+1$ .

Ainda, o processo $p=n_{p}-1$ deve computar

u_{j}=u^{0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j-1}-2\frac{h_{t}}{h_{x}^{2}}u^{% 0}_{j}+\frac{h_{t}}{h_{x}^{2}}u^{0}_{j+1},

(3.30)

para $j=i_{n_{p}-1}+1,\dotsc,f_{n_{p}-1}-1$ , com as condições de contorno

	$\displaystyle u^{0}_{i_{n_{p}-1}}=u^{0}_{f_{n_{p}-2}}$		(3.31)
	$\displaystyle u^{0}_{f_{n_{p}-1}}=0.$		(3.32)

O que requer a comunicação entre os processos $n_{p}-2$ e o $n_{p}-1$ .

Código: calor.cc

⬇

1#include <stdio.h>

2#include <math.h>

4// API MPI

5#include <mpi.h>

7int main (int argc, char** argv) {

8 // Inicializa o MPI

9 MPI_Init(NULL, NULL);

11 // numero total de processos

12 int world_size;

13 MPI_Comm_size(MPI_COMM_WORLD, &world_size);

15 // ID (rank) do processo

16 int world_rank;

17 MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

19 // parametros

20 size_t M = 1000;

21 size_t I = 10;

23 // tamanho dos passos discretos

24 double ht = 1e-3;

25 double hx = 1.0/I;

26 double cfl = ht/(hx*hx);

28 // malha espacial local

29 size_t ip = world_rank * int (I/world_size);

30 size_t fp = (world_rank+1) * int (I/world_size);

31 if (world_rank == world_size-1)

32 fp = I;

33 size_t my_I = fp - ip;

35 double x[my_I+1];

36 for (size_t j=0; j<=my_I; j++)

37 x[j] = (ip+j) * hx;

39 // solucao local

40 double u0[my_I+1], u[my_I+1];

42 // condicao inicial

43 for (size_t j=0; j<=my_I; j++) {

44 u0[j] = sin (M_PI * x[j]);

45 }

47 // condicoes de contorno de Dirichlet

48 if (world_rank == 0)

49 u[0] = 0.0;

50 if (world_rank == world_size-1)

51 u[my_I] = 0.0;

53 // auxiliares

54 double u00, u0I;

57 // iteracoes no tempo

58 for (size_t m=0; m<M; m++) {

60 if (world_rank == 0) {

61 MPI_Send (&u0[my_I], 1, MPI_DOUBLE,

62 1, 0, MPI_COMM_WORLD);

63 MPI_Recv (&u0I, 1, MPI_DOUBLE,

64 1, 0, MPI_COMM_WORLD,

65 MPI_STATUS_IGNORE);

66 } else if (world_rank < world_size-1) {

67 MPI_Recv (&u00, 1, MPI_DOUBLE,

68 world_rank-1, 0, MPI_COMM_WORLD,

69 MPI_STATUS_IGNORE);

70 MPI_Send (&u0[my_I], 1, MPI_DOUBLE,

71 world_rank+1, 0, MPI_COMM_WORLD);

73 MPI_Recv (&u0I, 1, MPI_DOUBLE,

74 world_rank+1, 0, MPI_COMM_WORLD,

75 MPI_STATUS_IGNORE);

76 MPI_Send (&u0[1], 1, MPI_DOUBLE,

77 world_rank-1, 0, MPI_COMM_WORLD);

78 } else {

79 MPI_Recv (&u00, 1, MPI_DOUBLE,

80 world_size-2, 0, MPI_COMM_WORLD,

81 MPI_STATUS_IGNORE);

82 MPI_Send (&u0[1], 1, MPI_DOUBLE,

83 world_size-2, 0, MPI_COMM_WORLD);

84 }

86 // atualizacao

87 u[1] = u0[1]

88 + cfl * u00

89 - 2*cfl * u0[1]

90 + cfl * u0[2];

91 for (size_t j=2; j<my_I; j++)

92 u[j] = u0[j]

93 + cfl * u0[j-1]

94 - 2*cfl * u0[j]

95 + cfl * u0[j+1];

96 if (world_rank < world_size-1)

97 u[my_I] = u0[my_I]

98 + cfl * u0[my_I-1]

99 - 2*cfl * u0[my_I]

100 + cfl * u0I;

101

102 // prepara nova iteracao

103 for (size_t j=0; j<=my_I; j++)

104 u0[j] = u[j];

105

106 }

107

108

109 // Finaliza o MPI

110 MPI_Finalize ();

111

112 return 0;

113}

3.5.1 Exercícios

Em revisão

E. 3.5.1.

O problema (3.15)-(3.17) tem solução analítica

u(t,x)=e^{-\pi^{2}t}\operatorname{sen}(\pi x).

(3.33)

E. 3.5.2.

Modifique o código calor.cc acima, de forma que a cada iteração $m$ no tempo, seja impresso na tela o valor do funcional

q(t)=\int_{0}^{1}u(t,x)\,dx.

(3.34)

Valide os resultados por comparação com a solução analítica (veja o Exercício (3.5.1)).

E. 3.5.3.

E. 3.5.4.

Modifique o código calor.cc acima, de forma que as comunicações sejam feitas de forma assíncrona, i.e. usando as rotinas MPI_Isend e MPI_Irecv. Há vantagens em se fazer esta modificação?

E. 3.5.5.

Faça um código MPI para computar uma aproximação da solução de

u_{t}-u_{xx}=0,\quad 0<x<1,t>0,

(3.35)

com condições de contorno periódicas

u(0,t)=u(1,t),t>0,

(3.36)

e condições iniciais

u(x,0)=\operatorname{sen}(\pi x),0\leq x\leq 1,

(3.37)

onde, $u=u(t,x)$ .

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Pedro H A Konzen

| | | |