Chương 5: Phương trình Poisson (Phần 2)

Xem lại Phần 1 • Trở về Mục lục cuốn sách

Phép biến đổi Fourier nhanh

Phương pháp được phác họa trong Mục [ĐK biên Dirichlet ở Phần 1] để giải phương trình Poisson trong 2 chiều với điều kiện biên Dirichlet đơn giản theo phương y đòi hỏi ta phải thực hiện rất nhiều lượt biến đổi Fourier cho hàm sin:

FjS = 2 / J ∑ k = 1J − 1fk sin(jkπ / J)  (ffta)

với j = 0, J, và các biến đổi Fourier ngược cho hàm sin:

fj = ∑ k = 1J − 1FkS sin(jkπ / J).   (fftb)

Ở đây, fj là giá trị của f(y) tại yj = jL / J. Như vậy, PT ([ffta]) giống như ([ffta1]) và ([ffta2]), còn PT ([fftb]) có thể được dùng để tái lập u(xi, yj) từ các Ui, j. Tương tự, phương pháp được vạch ra trong Mục [ĐK biên Neumann ở Phần 1] để giải hệ phương trình Poisson trong 2 chiều với các điều kiện biên Neumann đơn giản theo phương y đòi hỏi chúng ta phải thực hiện rất nhiều lượt biến đổi Fourier cho hàm cosin:

\displaystyle{ F_j^C = \frac{f_0}{J} + \frac{2}{J}\sum_{k=1}^{J-1} f_k\,\cos(j\,k\,\pi/J)  + \frac{(-1)^j\,f_J}{J} }

với j = 0, J, và các biến đổi Fourier ngược cho hàm cosin:

fj = ∑ k = 0JFkC cos(jkπ / J). 

Thật không may là để thực hiện trực tiếp các biến đổi như vậy, cần đến O(J2) phép tính số học, nghĩa là chúng cực kỳ tiêu tốn tài nguyên CPU. Tuy nhiên, có một thuật toán khéo léo dành riêng cho việc biến đổi Fourier mà chỉ tốn O(J lnJ) phép tính số học [vốn nhỏ hơn O(J2) đáng kể khi J lớn]. Thuật toán này có tên là biến đổi Fourier nhanh hay FFT.1

Chi tiết của thuật toán FFT vượt ra ngoài phạm vi khóa học này. Nói nôm na, thuật toán thực hiện việc chuyển đổi theo các giai đoạn nối tiếp dùng 2, 4, 8, 16, v.v. điểm nút. Ở khóa học này, ta sẽ dùng thư viện FFT thông dụng nhất có thể kiếm được trên mạng, thư viện fftw,2 để thực hiện tất cả các phép biến đổi Fourier cho hàm sin và cosin. Không may là thư viện fftw không trực tiếp thực hiện tính toán biến đổi Fourier các hàm sin và cosin.3 Thay vì vậy, nó tính biến đổi Fourier phức:

Fj = 1 / 2J ∑ k = 02J − 1fk exp( − i jkπ / J)

với j = 0, 2J − 1, và các biến đổi Fourier ngược phức:

fj = ∑ k = 02J − 1Fk exp( i jkπ / J). 

Lưu ý rằng fjFj có tính tuần hoàn theo j với chu kì 2 J. Cũng lưu ý rằng bộ số liệu gắn liền với phép biến đổi Fourier phức bao gồm số phần tử nhiều gấp hai lần số phần tử trong biến đổi sin hoặc cosin. Tuy nhiên, ta có thể dễ dàng chuyển từ phép biến đổi sin hoặc cosin sang biến đổi phức bằng cách kéo dài bộ số liệu. Vì vậy, đối với biến đổi sin ta viết:

\begin{array}{ccc}  f_{2J-j} &=& - f_j,\\[0.5ex]  F_{2J-j} &=& -F_j,  \end{array}

với j = 1, J − 1, trong trường hợp này

FjS = 2 i Fj

Tương tự, với biến đổi cosin ta viết:

\begin{array}{ccc}  f_{2J-j} &=& f_j,\\[0.5ex]  F_{2J-j} &=&F_j,  \end{array}

với j = 1, J − 1, trong trường hợp này

FjC = 2 Fj. 

Dưới đây là một loạt các đoạn chương trình để bọc các hàm trong thư viện fftw để thực hiện biến đổi Fourier cho các hàm sin và cosin.

// FFT.cpp 

// Set of functions to calculate Fourier-cosine and -sine transforms 
// of real data using fftw Fast-Fourier-Transform library. 
// Input/ouput arrays are assumed to be of extent J+1. 
// Uses version 2 of fftw library (incompatible with vs 3). 

#include <fftw.h> 
#include <blitz/array.h> 

using namespace blitz; 

// Calculates Fourier-cosine transform of array f in array F 
void fft_forward_cos (Array<double,1> f, Array<double,1>& F) 
{ 
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data 
  c_re (ff[0]) = f(0); c_im (ff[0]) = 0.; 
  c_re (ff[J]) = f(J); c_im (ff[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.; 
      c_re (ff[2*J-j]) = f(j); c_im (ff[2*J-j]) = 0.; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_FORWARD, FFTW_ESTIMATE); 
  fftw_one (p, ff, FF); 
  fftw_destroy_plan (p);  

  // Unload data 
  F(0) = c_re (FF[0]); F(J) = c_re (FF[J]);  
  for (int j = 1; j < J; j++) 
    { 
      F(j) = 2. * c_re (FF[j]); 
    } 

  // Normalize data 
  F /= 2. * double (J); 
} 

// Calculates inverse Fourier-cosine transform of array F in array f 
void fft_backward_cos (Array<double,1> F, Array<double,1>& f) 
{     
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data 
  c_re (FF[0]) = F(0); c_im (FF[0]) = 0.; 
  c_re (FF[J]) = F(J); c_im (FF[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (FF[j]) = F(j) / 2.; c_im (FF[j]) = 0.; 
      FF[2*J-j] = FF[j]; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_BACKWARD, FFTW_ESTIMATE); 
  fftw_one (p, FF, ff); 
  fftw_destroy_plan (p);  

  // Unload data  
  f(0) = c_re (ff[0]); f(J) = c_re (ff[J]);  
  for (int j = 1; j < J; j++) 
    { 
      f(j) = c_re (ff[j]); 
    } 
} 

// Calculates Fourier-sine transform of array f in array F 
void fft_forward_sin (Array<double,1> f, Array<double,1>& F) 
{   
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data  
  c_re (ff[0]) = 0.; c_im (ff[0]) = 0.;  
  c_re (ff[J]) = 0.; c_im (ff[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.; 
      c_re (ff[2*J-j]) = - f(j); c_im (ff[2*J-j]) = 0.; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_FORWARD, FFTW_ESTIMATE); 
  fftw_one (p, ff, FF); 
  fftw_destroy_plan (p);  

  // Unload data 
  F(0) = 0.; F(J) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      F(j) = - 2. * c_im (FF[j]); 
    }   

  // Normalize data 
  F /= 2. * double (J); 
} 

// Calculates inverse Fourier-sine transform of array F in array f 
void fft_backward_sin (Array<double,1> F, Array<double,1>& f) 
{  
  // Find J. Declare local arrays. 
  int J = f.extent(0) - 1; 
  int N = 2 * J; 
  fftw_complex ff[N], FF[N]; 

  // Load and extend data 
  c_re (FF[0]) = 0.; c_im (FF[0]) = 0.;  
  c_re (FF[J]) = 0.; c_im (FF[J]) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      c_re (FF[j]) = 0.; c_im (FF[j]) = - F(j) / 2.; 
      c_re (FF[2*J-j]) = 0.; c_im (FF[2*J-j]) = F(j) / 2.; 
    } 

  // Call fftw routine 
  fftw_plan p = fftw_create_plan (N, FFTW_BACKWARD, FFTW_ESTIMATE); 
  fftw_one (p, FF, ff); 
  fftw_destroy_plan (p);  

  // Unload data 
  f(0) = 0.; f(J) = 0.; 
  for (int j = 1; j < J; j++) 
    { 
      f(j) = c_re (ff[j]); 
    } 
}

Ví dụ tính toán 2 chiều cho bài toán Poisson 2 chiều

Dưới đây là một ví dụ đoạn chương trình giải bài toán Poisson 2 chiều có dùng các đoạn chương trình làm nhiệm vụ nghịch đảo ma trận ba đường chéo và bao bọc các hàm FFT nêu trên, cũng như thư viện Blitz++.

// Poisson2d.cpp 

// Function to solve Poisson's equation in 2-d: 

//  d^2 u / dx^2 + d^2 u / dy^2 = v  for  xl <= x <= xh  and  0 <= y <= L 

//  alphaL u + betaL du/dx = gammaL(y)  at x=xl 

//  alphaH u + betaH du/dx = gammaH(y)  at x=xh 

// In y-direction, either simple Dirichlet boundary conditions: 

//  u(x,0) = u(x,L) = 0 

// or simple Neumann boundary conditions: 

//  du/dy(x,0) = du/dy(x,L) = 0 

// Matrices u and v assumed to be of extent N+2, J+1. 
// Arrays gammaL, gammaH assumed to be of extent J+1. 

// Now, (i,j)th elements of matrices correspond to 

//  x_i = xl + i * dx    i=0,N+1 

//  y_j = j * L / J      j=0,J 

// Here, dx = (xh - xl) / (N+1) is grid spacing in x-direction. 

// Now, kappa = pi * dx / L 

// Finally, Neumann=0/1 selects Dirichlet/Neumann bcs in y-direction. 

#include <blitz/array.h> 

using namespace blitz; 

void fft_forward_cos (Array<double,1> f, Array<double,1>& F); 
void fft_backward_cos (Array<double,1> F, Array<double,1>& f); 
void fft_forward_sin (Array<double,1> f, Array<double,1>& F); 
void fft_backward_sin (Array<double,1> F, Array<double,1>& f); 
void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c,  
                  Array<double,1> w, Array<double,1>& u); 

void Poisson2D (Array<double,2>& u, Array<double,2> v,  
                double alphaL, double betaL, Array<double,1> gammaL,  
                double alphaH, double betaH, Array<double,1> gammaH, 
                double dx, double kappa, int Neumann) 
{ 
  // Find N and J. Declare local arrays. 
  int N = u.extent(0) - 2; 
  int J = u.extent(1) - 1; 
  Array<double,2> V(N+2, J+1), U(N+2, J+1); 
  Array<double,1> GammaL(J+1), GammaH(J+1); 

  // Fourier transform boundary conditions 
  if (Neumann) 
    { 
      fft_forward_cos (gammaL, GammaL); 
      fft_forward_cos (gammaH, GammaH); 
    } 
  else 
    { 
      fft_forward_sin (gammaL, GammaL); 
      fft_forward_sin (gammaH, GammaH); 
    } 

  // Fourier transform source term 
  for (int i = 1; i <= N; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      for (int j = 0; j <= J; j++) In(j) = v(i, j); 

      if (Neumann) 
        fft_forward_cos (In, Out); 
      else 
        fft_forward_sin (In, Out); 

      for (int j = 0; j <= J; j++) V(i, j) = Out(j); 
    } 

  // Solve tridiagonal matrix equations 
  if (Neumann) 
    { 
      for (int j = 0; j <= J; j++) 
        {   
          Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2), uu(N+2); 

          // Initialize tridiagonal matrix 
          for (int i = 2; i <= N; i++) a(i) = 1.; 
          for (int i = 1; i <= N; i++) 
            b(i) = -2. - double (j * j) * kappa * kappa; 
          b(1) -= betaL / (alphaL * dx - betaL); 
          b(N) += betaH / (alphaH * dx + betaH); 
          for (int i = 1; i <= N-1; i++) c(i) = 1.; 

          // Initialize right-hand side vector 
          for (int i = 1; i <= N; i++) 
            w(i) = V(i, j) * dx * dx; 
          w(1) -= GammaL(j) * dx / (alphaL * dx - betaL); 
          w(N) -= GammaH(j) * dx / (alphaH * dx + betaH); 

          // Invert tridiagonal matrix equation 
          Tridiagonal (a, b, c, w, uu); 
          for (int i = 1; i <= N; i++) U(i, j) = uu(i); 
        } 
    } 
  else 
    { 
      for (int j = 1; j < J; j++) 
        {  
          Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2), uu(N+2); 

          // Initialize tridiagonal matrix 
          for (int i = 2; i <= N; i++) a(i) = 1.; 
          for (int i = 1; i <= N; i++) 
            b(i) = -2. - double (j * j) * kappa * kappa; 
          b(1) -= betaL / (alphaL * dx - betaL); 
          b(N) += betaH / (alphaH * dx + betaH); 
          for (int i = 1; i <= N-1; i++) c(i) = 1.; 

          // Initialize right-hand side vector 
          for (int i = 1; i <= N; i++) 
            w(i) = V(i, j) * dx * dx; 
          w(1) -= GammaL(j) * dx / (alphaL * dx - betaL); 
          w(N) -= GammaH(j) * dx / (alphaH * dx + betaH); 

          // Invert tridiagonal matrix equation 
          Tridiagonal (a, b, c, w, uu); 
          for (int i = 1; i <= N; i++) U(i, j) = uu(i); 
        } 

      for (int i = 1; i <= N  ; i++)  
        { 
          U(i, 0) = 0.; U(i, J) = 0.; 
        } 
    } 

  // Reconstruct solution via inverse Fourier transform 
  for (int i = 1; i <= N; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      for (int j = 0; j <= J; j++) In(j) = U(i, j); 

      if (Neumann) 
        fft_backward_cos (In, Out); 
      else 
        fft_backward_sin (In, Out); 

      for (int j = 0; j <= J; j++) u(i, j) = Out(j); 
    } 

  // Calculate i=0 and i=N+1 values   
  for (int j = 0; j <= J; j++) 
    { 
      u(0, j) = (gammaL(j) * dx - betaL * u(1, j)) / 
        (alphaL * dx - betaL); 
      u(N+1, j) = (gammaH(j) * dx + betaH * u(N, j)) / 
        (alphaH * dx + betaH); 
    } 
}

Một ví dụ giải phương trình Poisson 2 chiều

Bây giờ ta hãy dùng kĩ thuật nêu trên để giải phương trình Poisson trong không gian hai chiều. Giả sử rằng số hạng nguồn là

v(x, y) = 6 xy (1 − y) − 2 x3

với 0 ≤ x ≤ 10 ≤ y ≤ 1. Các điều kiện biên tại x = 0αl = 1, βl = 0, và γl = 0 [xem PT ([e527a])], trong đó các điều kiện biên tại x = 1αh = 1, βh = 0, và γh = y (1 − y) [xem PT ([e528a])]. Các điều kiện biên Dirichlet đơn giản u(x, 0) = u(x, 1) = 0 được áp dụng tại y = 0y = 1. Dĩ nhiên là bài toán này có thể giải được theo cách giải tích để cho ta

u(x, y) = y (1 − y) x3. 

Các Hình [p2da] và [p2db] so sánh giữa các nghiệm giải tích và sai phân hữu hạn cho N = J = 64. Có thể thấy rằng nghiệm sai phân hữu hạn khá trùng khớp với nghiệm giải tích.

p2daHình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trục x tại y = 0.5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2da]

p2db

Hình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trục y tại x = 0.5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2db]

Xét ví dụ thứ hai, giả sử rằng số hạng nguồn là

v(x, y) =  − 2 (2 y3 − 3 y2 + 1) + 6 (1 − x2) (2 y − 1)

với 0 ≤ x ≤ 10 ≤ y ≤ 1. Các điều kiện biên tại x = 0αl = 1, βl = 0, và γl = 2 y3 − 3 y2 + 1 [xem PT ([e527a])], còn các điều kiện biên tại x = 1αh = 1, βh = 0, và γh = 0 [xem PT ([e528a])]. Các điều kiện biên Neumann đơn giản là u(x, 0) / ∂y = ∂u(x, 1) / ∂y = 0 được áp dụng tại y = 0y = 1. Dĩ nhiên là bài toán này có thể giải được theo cách giải tích để cho ta

u(x, y) = (1 − x2) (2 y3 − 3 y2 + 1). 

Các Hình [p2dc] và [p2dd] cho thấy sự so sánh giữa các nghiệm giải tích và sai phân hữu hạn cho N = J = 64. Có thể thấy rằng nghiệm sai phân hữu hạn khá trùng khớp với nghiệm giải tích.

p2dc

Hình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trục x tại y = 0. 5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2dc]

p2dd

Hình. Nghiệm của phương trình Poisson trong không gian hai chiều với các điều kiện biên Dirichlet đơn giản theo phương y. Nghiệm được vẽ theo trục y tại x = 0.5. Đường chấm chấm (bị che khuất) biểu thị nghiệm giải tích, còn các hình tam giác biểu thị nghiệm sai phân hữu hạn cho N = J = 64. [p2dd]

Ví dụ tính toán trường tĩnh điện 2 chiều

Ta hãy làm một ví dụ tính toán trường tĩnh điện 2 chiều. Xét một dây dẫn mang điện chạy song song trục của một đường kênh dẫn rỗng hình hộp chữ nhật. Giả sử rằng các đỉnh của hình hộp này nằm ở (x, y) = (0, 0), (0, 1), (1, 0), và (1, 1). Cũng giả sử rằng dây dẫn mang điện tích đều bằng 1 đơn vị trên mỗi đơn vị dài của dây. Điện thế φ(x, y) bên trong kênh dẫn thỏa mãn [xem PT ([estat])]

2φ(x, y) / ∂x2 + ∂2φ(x, y) / ∂y2 = v(x, y) =  − δ(x − x0) δ(y − y0),   (estat1)

troing đó (x0, y0) là các tọa độ của dây. Ở đây để cho tiện ta đã chuẩn hóa các đơn vị sao cho hệ số ε0 được hấp thụ trong quá trình chuẩn hóa. Giả sử rằng hộp được nối đất, điện thế sẽ tuân theo các điều kiện biên Dirichlet φ = 0 tại x = 0, x = 1, y = 0, và y = 1. Ta cần tìm nghiệm trong vùng 0 ≤ x ≤ 10 ≤ y ≤ 1.

coulomb1

Hình. Đồ thị đường đồng mức của điện thế hình thành do một dây dẫn được tích điện đặt ở tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại (x, y) = (0. 5, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub1]

Lưu ý rằng khi làm rời rạc PT ([estat1]) vế phải trở thành

v(xi, yj) =  − 1 / (δxδy)

tại điểm nút lưới sát dây dẫn nhất, với v(xi, yj) = 0 trên các điểm nút lưới còn lại. Ở đây, δxδy lần lượt là các khoảng cách nút lưới theo các phương xy.

coulomb2

Hình. Đồ thị véc-tơ thể hiện hướng của điện trường hình thành do một dây dẫn được tích điện đặt ở tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại (x, y) = (0. 5, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub2]

Các Hình [coub1] và [coub2] biểu diễn điện thế φ(x, y) và điện trường E =  − ∇φ hình thành do một dây dẫn được tích điện đặt ở tâm một kênh dẫn chữ nhật, tức là tại (x0, y0) = (0. 5, 0. 5). Việc tính toán được thực hiện với chương trình giải bài toán Poisson 2 chiều nêu trên với N = J = 128.

coulomb3

Hình. Đồ thị đường đồng mức của điện thế hình thành do một dây dẫn được tích điện đặt lệch tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại (x, y) = (0. 25, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub3]

Các Hình [coub3] và [coub4] biểu diễn điện thế φ(x, y) và điện trường E =  − ∇φ hình thành do một dây dẫn được tích điện đặt lệch khỏi tâm kênh dẫn chữ nhật, tức là tại (x0, y0) = (0. 25, 0. 5). Việc tính toán được thực hiện với chương trình giải bài toán Poisson 2 chiều nêu trên với N = J = 128.

coulomb4

Hình. Đồ thị véc-tơ thể hiện hướng của điện trường hình thành do một dây dẫn được tích điện đặt lệch tâm một kênh dẫn chữ nhật được nối đất. Dây được đặt tại (x, y) = (0. 25, 0. 5), còn các vách kênh dẫn tại x = 0, x = 1, y = 0, và y = 1. Tính toán được thực hiện với N = J = 128. [coub4]

Bài toán 3 chiều

Những kĩ thuật được thảo luận trong các Mục [fftsine] và [fftcos] để giải phương trình Poisson trong không gian hai chiều với một lớp các điều kiện biên hạn hẹp có thể được dễ dàng mở rộng ra cho không gian 3 chiều. Trong trường hợp 3 chiều, ta cần phải thực hiện biến đổi Fourier theo hai chiều (chẳng hạn theo các phương yz) để giản hóa bài toán về một ma trận ba đường chéo nhưng các phương trình không ràng buộc nhau. Những phương trình này có thể được nghịch dảo theo cách thông thường, và nghiệm được tái lập lại bằng biến đổi Fourier ngược hai lần.


  1. Xem Numerical recipes in C: the art of scientific computing, W.H. Press, S.A. Teukolsky, W.T. Vettering, and B.R. Flannery (Cambridge University Press, Cambridge, England, 1992), Mục 12.2.
  2. Xem http://www.fftw.org
  3. Điều này đúng cho thư viện phiên bản 2 (dùng cho khóa học này), nhưng phiên bản 3 thì khác.
Advertisements

2 phản hồi

Filed under Vật lý tính toán

2 responses to “Chương 5: Phương trình Poisson (Phần 2)

  1. Pingback: Chương 1: Giới thiệu chung | Blog của Chiến

  2. Pingback: Tìm hiểu về biến đổi Fourier FFT/IFFT và phương pháp tính toán | Vi-et Spaces

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Log Out / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Log Out / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Log Out / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Log Out / Thay đổi )

Connecting to %s