периодические граничные условия - конечные разности

Привет, у меня есть код ниже, который решает нелинейно связанные PDE. Однако мне нужно реализовать периодические граничные условия. Периодические граничные условия беспокоят меня, что я должен добавить в свой код, чтобы обеспечить соблюдение периодических граничных условий? Обновлено на основе приведенных ниже модульных арифметических предложений.

Обратите внимание, t>=0 и x находится в интервале [0,1]. Вот связанные уравнения, ниже я привожу свой код

уравнения

где а, b > 0.

Это начальные условия, но теперь мне нужно наложить периодические граничные условия. Математически их можно записать как u(0,t)=u(1,t) и du(0,t)/dx=du(1,t)/dx, то же самое верно для f(x,t). du/dx, которые у меня есть для граничных условий, на самом деле должны быть частными производными.

Мой код ниже

program coupledPDE 

integer, parameter :: n = 10, A = 20 
real, parameter :: h = 0.1, k = 0.05 
real, dimension(0:n-1) :: u,v,w,f,g,d 
integer:: i,m 
real:: t, R, x,c1,c2,c3,aa,b 

R=(k/h)**2.
aa=2.0
b=1.0
c1=(2.+aa*k**2.-2.0*R)/(1+k/2.)
c2=R/(1.+k/2.)
c3=(1.0-k/2.)/(1.0+k/2.)
c4=b*k**2./(1+k/2.)


do i = 0,n-1 !loop over all space points except 0 and n
  x = real(i)*h    
  w(i) = z(x)  !u(x,0)
  d(i) = z(x)  !f(x,0)
end do


do i=0,n-1
  ip=mod(i+1,n)
  il=modulo(i-1,n)
  v(i) = (c1/(1.+c3))*w(i) + (c2/(1.+c3))*(w(ip)+w(il)) -(c4/(1.+c3))*w(i)*((w(i))**2.+(d(i))**2.)    !\partial_t u(x,0)=0
  g(i) = (c1/(1.+c3))*d(i) + (c2/(1.+c3))*(d(ip)+d(il)) -(c4/(1.+c3))*d(i)*((w(i))**2.+(d(i))**2.)    !\partial_t f(x,0)=0
end do

do m=1,A 

   do i=0,n-1
       ip=mod(i+1,n)
       il=modulo(i-1,n)
       u(i)=c1*v(i)+c2*(v(ip)+v(il))-c3*w(i)-c4*v(i)*((v(i))**2.+(g(i))**2.) 
       f(i)=c1*g(i)+c2*(g(ip)+g(il))-c3*d(i)-c4*g(i)*((v(i))**2.+(g(i))**2.) 
   end do 
     print*, "the values of u(x,t+k) for all m=",m
   print "(//3x,i5,//(3(3x,e22.14)))",m,u   

  do i=0,n-1
   w(i)=v(i)
   v(i)=u(i)
   d(i)=g(i)
   t=real(m)*k
   x=real(i)*k
  end do

end do


end program coupledPDE

function z(x)
real, intent(in) :: x
real :: pi
pi=4.0*atan(1.0)
z = sin(pi*x)
end function z

Спасибо за чтение, если я должен переформатировать свой вопрос более правильным образом, пожалуйста, дайте мне знать.


fortran math numerical-methods pde fortran90
person Jeff Faraci    schedule 09.03.2016    source источник
comment
Просто скопируйте значения от n до 0 и от 1 до n+1 на каждом временном шаге. Это все.   -  person Vladimir F    schedule 10.03.2016
comment
Конечно, чище объявить все только на сетке 0,h,2h,...,1-h, а затем использовать простую арифметику по модулю для сдвига индексов. О периодических границах заботятся автоматически. Еще лучше было бы использовать целые массивы и функцию CSHIFT.   -  person RussF    schedule 10.03.2016
comment
@VladimirF Спасибо за помощь, однако я не понимаю, что вы имеете в виду. Где я копирую эти значения и какие значения вы имеете в виду? Должен ли я добавить код, когда я запускаю свою временную петлю, или? (Я понимаю, что это действительно простая идея для реализации, извините за путаницу с моей стороны). Спасибо еще раз!   -  person Jeff Faraci    schedule 10.03.2016
comment
@RussF Спасибо и за ответ. Как объявить это на сетке? У меня реальная проблема с пониманием того, как это сделать. Спасибо за помощь   -  person Jeff Faraci    schedule 10.03.2016
comment
Просто объявите все свои массивы как v[0:n-1] и т. д. В ваших циклах вычислите ip=mod(i+1,n) и замените ссылки на i+1 на ip. Аналогично для i-1. Если вы определяете свои циклы как do i = 0,n-1, все обрабатывается естественным образом. Если вам действительно нужно значение x=1, просто скопируйте значение x=0 после цикла i.   -  person RussF    schedule 10.03.2016
comment
О, да, у вас ошибка границ в первом цикле. Это должно быть do i=0,n, а не do i=0,n+1   -  person RussF    schedule 10.03.2016
comment
@RussF Я исправил ошибку границ в первом цикле. Спасибо. Я изменил размеры всех массивов, как вы предлагаете. Я просто хочу уточнить то, что вы сказали, я использовал модульную арифметику для двух циклов дискретизации для i + 1, i-1. Также все мои циклы do, кроме первого, были изменены с i=0 на n-1. Мой вопрос заключается в следующем: почему нам не нужна модульная арифметика для факторов v(i), g(i) и т. д. Еще раз спасибо! Кроме того, когда я сделал это, в моей команде записи я заметил, что теперь uf не равен нулю.   -  person Jeff Faraci    schedule 10.03.2016
comment
Покажи, что ты сделал. Приведенный выше код не будет компилироваться из-за ошибок. Также убедитесь, что вы используете функцию modulo для i-1, а не mod. Вы не обновляете v, w, g или d, поэтому ничего не изменится.   -  person RussF    schedule 10.03.2016
comment
@RussF Спасибо, исправлено. Я только что обновил свой код в сообщении выше, чтобы показать вам, что я сделал, используя mod, modulo и изменяя циклы do. Я пропустил часть кода, в которой изначально обновлял v,w,g,d; но теперь это в коде выше. Теперь я замечаю, что на каждом временном шаге, если я печатаю решение, 2-я и 10-я точки совпадают, поэтому я вижу периодичность. Однако правильно ли я это реализовал? Результаты отличаются от того, когда я использовал метод призрака, как описано ниже. Я предпочитаю использовать модульную арифметику, так как мне нравится делать все на сетке, как вы предложили. Спасибо!   -  person Jeff Faraci    schedule 10.03.2016

Ответы (1)


arrow_upward
2
arrow_downward

Одним из вариантов граничных условий при дискретизации PDE является использование призрачных (гало) ячеек (узлов сетки). Возможно, он не самый умный для периодического BC, но его можно использовать для всех других типов граничных условий.

Итак, вы объявляете свои массивы как

real, dimension(-1:n) :: u,v,w,f,g,d

но вы решаете свое УЧП только в точках 0..n-1 (точка n совпадает с точкой 0). Вы также можете сделать из 1..n и объявить массивы из 0..n+1.

Затем вы устанавливаете

 u(-1) = u(n-1)

а также

 u(n) = u(0)

и то же самое для всех остальных массивов.

На каждом временном шаге вы снова устанавливаете это для u и f или всех других полей, которые изменяются во время решения:

do m=1,A 
   u(-1) = u(n-1)
   u(n) = u(0)
   f(-1) = f(n-1)
   f(n) = f(0)

   do i=0,n-1 !Discretization equation for all times after the 1st step
       u(i)=...
       f(i)=...
   end do 
end do

Все вышеизложенное предполагало явную временную дискретизацию и пространственную дискретизацию с конечными разностями и предполагало, что x(0) = 0 и x(n) = 1 являются вашими граничными точками.

person Vladimir F    schedule 10.03.2016
comment
Большое спасибо за помощь! Прежде чем я смогу комментировать дальше, мне нужно пройти через все это и подумать об этом, а также реализовать это. Я вернусь позже. - person Jeff Faraci; 10.03.2016
comment
Я нарисовал сетку и теперь понимаю, что вы имеете в виду, почему u(-1)=u(n-1) и u(n)=u(0), поэтому я понимаю, почему вы изменили размеры массива. Это очень полезно! Для двух циклов дискретизации я применил эти условия ко всем переменным u,f,v,g. Наконец, когда я распечатываю свое решение на каждом временном шаге, теперь я вижу, что действительно 1-е и 11-е значения совпадают! Спасибо за подробное объяснение, пожалуйста, дайте мне знать, если я сделал что-то неправильно. У меня есть один вопрос, этот метод все еще стабилен/рекомендуется в более высоких пространственных измерениях? Спасибо! - person Jeff Faraci; 10.03.2016
comment
@Integrals Я лично использую этот метод в вычислениях 3D CFD. Преимущество в том, что ваши вычислительные операторы проще и везде одинаковы. Недостатком является то, что вам нужно часто копировать границы, и это немного сложнее для неявных методов. Возможен и другой вариант с modulo. Я нахожу клетки-призраки более общими. - person Vladimir F; 10.03.2016
comment
Спасибо за информативный ответ, это очень полезно для меня. - person Jeff Faraci; 10.03.2016
comment
Привет, извините, что спрашиваю вас о чем-то еще здесь, если я должен опубликовать это как вопрос, пожалуйста, дайте мне знать. Мне было интересно, правильно ли мое обобщение вашего кода выше для двумерных периодических граничных условий: u(0,:)=u(nx,:); u(-1,:)=u(nx-1,:); и(:,0)=и(:,ни); u(:,-1)=u(:,ny-1) где nx,ny - это просто общее количество шагов по x и y, у меня квадратная сетка, поэтому они равны. Спасибо за вашу помощь. (Я бы наложил это на все массивы, а также добавил бы в каждый из двух циклов дискретизации). Я понимаю, почему вы сказали, что вам нужно часто копировать границы. - person Jeff Faraci; 12.03.2016