Объявление пространственных точек реалами вместо целых - PullRequest
0 голосов
/ 09 октября 2018

Я попытался продемонстрировать ниже в примере кода, что я пытаюсь сделать.У меня есть куча условий if / else, и я хочу правильно ввести их, когда я компилирую и запускаю код.Однако я не буду вводить те, где i, j - действительное число, поскольку i, j определяется как целое число в начале кода.Как я могу смоделировать эту проблему по-другому, чтобы я мог ввести те условия if / else?Спасибо

В частности, я не хочу, чтобы пространственные точки, i, j были целыми числами.Вместо этого я хочу, чтобы они были реальными значениями, которые соответствуют любой позиции в пространстве, например: (x, y) = (5.31, 5.31).

program sample

integer :: i,j !spatial points
real, parameter :: DeltaX = 0.1
integer, parameter :: n = 10
real, dimension(-n:n, -n:n) :: u,f

do j = -n,n
do i = -n,n

   f(i,j) = sin(i*DeltaX+j*DeltaX)

end do
end do

do i = -n+1,n-1 !this loop is only over integers
do j = -n+1,n-1 !this loop is only over integers

  if (i == 6 .AND. j == 6 .AND. i == j) then  

    PRINT*, 'i,j', i,j

    u(i,j) = (f(i+1,j+1)+f(i-1,j-1)-f(i,j))/DeltaX**2

  else if ( i == 8 .AND. j == 8. .AND. i == j) then 

    PRINT*, 'i,j ', i,j

    u(i,j) = (f(i+1,j+1)+f(i-1,j-1)-f(i,j))/DeltaX**2

  else if ( i == 5.31 .AND. j == 5.31 .AND. i == j) then 

    PRINT*, 'i,j = ', i,j !I will never enter this if statement though

    u(i,j) = (f(i+1,j+1)+f(i-1,j-1)-f(i,j))/DeltaX**2

  else if ( i == -6.87 .AND. j == -6.87 .AND. i == j) then 

    PRINT*, 'i,j = ', i,j !I will never enter this if statement either

    u(i,j) = (f(i+1,j+1)+f(i-1,j-1)-f(i,j))/DeltaX**2

  else 

   u(i,j) = (f(i+1,j)+f(i-1,j)+f(i,j+1)+f(i,j-1)-4*f(i,j))/DeltaX**2

  end if

end do
end do

end program

1 Ответ

0 голосов
/ 09 октября 2018

Похоже, что вы хотите решить дифференциальное уравнение в частных производных или какую-то другую математическую задачу, которая возникает путем дискретизации непрерывного уравнения в конечный набор точек.Точки часто находятся на ортогональной сетке в направлениях x и y , но сетка может быть повернута, наклонена, деформирована или даже полностью неструктурирована.Полезно помнить, что поскольку это подразумевает, что вам нужно что-то другое, что вы пробовали.

Для структурированных сеток вы можете использовать массивы Fortran для хранения значений переменных в каждой точке данных.

real ::  u(1:nx,1:ny), f(1:nx, 1:ny)

Вы можете сделать то же самое для x и y координат каждой точки.

Для простых ортогональных сеток в x и y вы можете просто использовать1D массивы для x и y координат

real :: x(1:nx), y(1:ny)

Смысл этого состоит в том, что каждая точка с данным x index i имеет значение x -координировать x(i) и аналогично для y(j).

Когда ваши вычисления зависят от координат, вы можете использовать их просто как

do j = 1, ny
  do i = 1, nx
    f(i,j) = func(x(i), y(j))
  end do
end do

И для ваших операторов if:

do j = 1, ny
  do i = 1, nx
    if (x(i)>something .and. x(i)<something_else) then
      u(i,j) = some_expression with u(i+-1,j+-1) and f(i+-1,j+-1)
    else ...
    end if
  end do
end do

Обратите внимание, что вы никогда не должны сравнивать значения с плавающей запятой (Fortran real) на равенство, подобное if (x==5.31).Числа с плавающей точкой неточны, и x может легко быть 5.3100000001 вместо 5.31.


Вы можете найти бесчисленное множество примеров в Интернете и даже здесь, в Переполнении стека.Просто поищите в Интернете «уравнение Фортрана Пуассона», и вы найдете огромное количество простых примеров кода для метода Якоби и других более сложных методов.

Помните, что для научных вычислений существует специальный сайт StackExchange,Мы могли бы сразу отправить вас туда, но ваша первая версия вопроса сильно отличалась и не очень подходила для этого сайта.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...