Fortran 90-95 코드에서 포화되지 않은 토양을 통한 물의 흐름을 시뮬레이트하는 데있어 잘못된 점은 무엇입니까?

Question

나는이 코드의 잘못된 점을 찾으려고 여러 해 동안 노력 해왔다. 포화되지 않은 토양을 통한 물의 흐름을 모델링하는 데 사용됩니다. 방정식 시스템은 삼중 항 행렬의 형태로 토마스 알고리즘으로 해결됩니다. 나는 해결책을 가지고 있으며, 코드는 그것을 표현하지 못한다. 예를 들어, 노드 A는 aprox -100 cm의 초기 조건에서 aprox -20 cm로가는 곡선이어야합니다. 그것은 긴 코드입니다,하지만 누군가가 저를 도왔다면 대단히 감사 할 것입니다.

program EcuacionRichards 

implicit none 

!Declaring variables 

integer, parameter :: nodos = 100 
integer :: i, it, max_it, nodo_a, nodo_b, nodo_c, nodo_d, it_bajo, it_alto 
double precision, dimension(1:nodos) :: H, H_ant, C, K, theta, theta_ant, aa, bb, cc, dd, rr, th_ant 
double precision :: dz, zbot, tfin, dt, rz, Ksup, Kinf, t, th_lisimetro, h_lisimetro 
double precision :: q_ent, tol_h, tol_th, cambio_h, cambio_th 
double precision :: mult_alto, mult_bajo, maxdt, mindt, qlibre 
logical lisimetro 

!Hydraulic Parameters 
double precision :: theta_sat=0.43  !cm/cm 
double precision :: theta_res=0.078  !cm/cm 
double precision :: alpha=0.0325  !1/cm 
double precision :: n=1.346 
double precision :: m 
double precision :: K_sat=86.4   !cm/d 

!Grid and iteration parameters 
lisimetro=.true. 
dt=0.01    !days 
zbot=160   !depth of the column in cm 
dz=zbot/nodos  !cm 
tfin=30    !days 
max_it=500   !max number of Picard iterations 
tol_h=0.1   !tolerance for H iteration, cm 
tol_th=0.001  !tolerance for theta iteration, 1/1 
it_bajo=3   !minimum recommended number of iterations 
it_alto=7   !maximum recommended number of iterations 
mult_bajo=1.3  !time multiplicator for low iterations 
mult_alto=0.7  !time multiplicator for low iterations 
maxdt=0.5   !max value for dt 
mindt=0.001   !min value for dt 
m=1-1/n 

!Initializing other variables 
th_lisimetro=0.32 
h_lisimetro=HfTH(th_lisimetro) 

nodo_a=nodos 
nodo_b=2*nodos/3 
nodo_c=nodos/3 
nodo_d=1 

!*********Initial Conditions************************************************************ 

call theta_ini(theta,nodos) !Fill array with initial moisture values 
do i=1,nodos 
    H(i)=HfTH(theta(i)) 
    call actualiza(H(i), theta(i), C(i), K(i)) 
end do 

!************* OPEN WRITING FILES ************************************************ 
open(unit=1,file='succion2.txt') 
open(unit=2,file='humedad2.txt') 
open(unit=3,file='conducti2.txt') 
open(unit=4,file='parametr2.txt') 
write(4,'("dt(días) =",f7.4)') dt 
write(4,'("dz(cm) =",f7.4)') dz 
write(4,'("nodos =",i5)') nodos 
write(4,'("altura(cm) =",f8.3)') zbot 
write(4,'("tfin(días) =",f7.2)') tfin 
write(4,'("theta_sat =",f7.4)') theta_sat 
write(4,'("theta_res =",f7.4)') theta_res 
write(4,'("K_saturada =",g11.3)') K_sat 
write(4,'("n =",f7.4)') n 
write(4,'("m =",f7.5)') m 
write(4,'("alpha =",f7.5)') alpha 
write(4,'("max_it =",i4)') max_it 
close(4) 
write(1,*) "T(días) H_a(cm) H_b(cm) H_c(cm) H_d(cm)" 
write(2,*) "T(días) th_a(cm) th_b(cm) th_c(cm) th_d(cm)" 
write(3,*) "T(días) K_a(cm/d) K_b(cm/d) K_c(cm/d) K_d(cm/d)" 


!*************TIME LOOP********************************************************************************************** 
t=0.d0 
do while ((t.le.tfin).and.(dt.gt.0)) 
    rz=dz/dt 
    t=t+dt 
    theta_ant=theta !Previous time 
    !Water flow that enters at the top (constant) 
    q_ent=0.1  !cm/dia 
!*************  PICARD LOOP    ****************************************** 
Picard:do it=1,max_it 

      if(it.eq.max_it) pause "MAXIMUM ITERATIONS REACHED" 

      !Interior Nodes 
      do i=2, nodos-1 
       Ksup=2*(K(i+1)*K(i))/(K(i+1)+K(i)) 
       Kinf=2*(K(i-1)*K(i))/(K(i-1)+K(i)) 
       aa(i)=-Kinf/dz !K(i-1/2) 
       cc(i)=-Ksup/dz !K(i+1/2) 
       bb(i)=rz*C(i)-aa(i)-cc(i) 
       rr(i)=rz*C(i)*h(i)-rz*(theta(i)-theta_ant(i))+Ksup-Kinf 
      end do 

      !Inferior Node 
      if (lisimetro) then 
       !Changing inferior node 
       if (theta(1).lt.th_lisimetro) then 
        !Water flow 0, Neumann 
        Ksup=2*(K(1)*K(2))/(K(1)+K(2)) 
        aa(1)=0 
        cc(1)=-Ksup/dz 
        bb(1)=-cc(1) 
        rr(1)=Ksup 
       else 
        !H(1)=0 condition, Dirichlet 
        Ksup=2*(K(1)*K(2))/(K(1)+K(2)) 
        aa(1)=0 
        bb(1)=1 
        cc(1)=0 
        rr(1)=h_lisimetro 
        aa(2)=0 
        rr(2)=rr(2)+Ksup/dz*(h_lisimetro) 
       end if 
      else 
       !Inferior node, free drainage, Neumann 
       Ksup=2*(K(1)*K(2))/(K(1)+K(2)) 
       qlibre=-K(1) 
       aa(1)=0 
       cc(1)=-Ksup/dz 
       bb(1)=-cc(1) 
       rr(1)=Ksup+qlibre 
      end if 

      !Superior node, known water flow 
      Kinf=2*(K(nodos)*K(nodos-1))/(K(nodos)+K(nodos-1)) 
      aa(nodos)=-Kinf/dz 
      cc(nodos)=0 
      bb(nodos)=0.5*rz*C(nodos)-aa(nodos) 
      rr(nodos)=0.5*rz*C(nodos)*h(nodos)-0.5*rz*(theta(nodos)-theta_ant(nodos))-Kinf-q_ent 

      call tridiag(aa,bb,cc,rr,dd,nodos) 

      !Suction modification and H functions actualization 
      h_ant=h 
      th_ant=theta !Save iteration 
      h=dd   !Advance to next iteration 
      do i=1,nodos 
       call actualiza(H(i),theta(i), C(i), K(i)) 
      end do 

      !End of iterations condition 
      cambio_h=maxval(dabs(h-h_ant)) 
      cambio_th=maxval(dabs(theta-th_ant)) 

      if((cambio_h.lt.tol_h).and.(cambio_th.lt.tol_th)) then 

       if(.true.) then !(t.eq.tprint) 
       write (1,'(f8.3,f9.3,f9.3,f9.3,f9.3)') t,H(nodo_a),H(nodo_b),H(nodo_c),H(nodo_d) 
       write (2,'(f8.3,f7.4,f7.4,f7.4,f7.4)') t,theta(nodo_a),theta(nodo_b),theta(nodo_c),theta(nodo_d) 
       write (3,'(f8.3,g11.4,g11.4,g11.4,g11.4)') t,k(nodo_a),k(nodo_b),k(nodo_c),k(nodo_d) 
       end if 

       if (it.lt.it_bajo) dt=min(dt*mult_bajo,maxdt) 
       if (it.gt.it_alto) dt=max(dt*mult_alto,mindt) 

       exit Picard 

      else 
       cycle Picard 
      end if 
     end do Picard !Picard loop end 
     if ((tfin-t).le.1E-4) t=huge(1.d0) 
end do 
!Time Loop End*************************************************************** 
!******** Close files 
close(1) 
close(2) 
close(3) 

!********END OF PROGRAM********************************************************** 
!****************************************************************************** 
!Subroutines and functions 
contains 

!Initial moistures assignment 
subroutine theta_ini(theta,nodos) 
integer :: nodos 
double precision, dimension(1:nodos) :: theta 
integer i 
do i=1, nodos 
    theta(i)=0.30 
end do 
end subroutine theta_ini 

!Subroutine that actualizes salues according to pressure 
subroutine actualiza(p,theta,c,k) 
    double precision p, theta, c, k 
    double precision se, te 
    if(p.lt.0) then   
        te=1+(-alpha*p)**n 
        se=te**(-m) 
        theta=theta_res+(theta_sat-theta_res)*se 
        K=K_sat*se**(0.5)*(1-(1-se**(1/m))**m)**2 
        c=((alpha**n)*(theta_sat-theta_res)*n*m*(-p)**(n-1))/(te**(m+1)) !d(theta)/dh 
    else 
        theta=theta_sat 
        K=K_sat 
        c=0 
    end if 
    return 
end subroutine actualiza 

!Tridiag(alpha,beta, gamma, Resto, delta, nodos) 
     subroutine tridiag(a,b,c,d,x,n) 
     implicit none 
!  a - sub-diagonal (means it is the diagonal below the main diagonal) 
!  b - the main diagonal 
!  c - sup-diagonal (means it is the diagonal above the main diagonal) 
!  d - right part 
!  x - the answer 
!  n - number of equations 

     integer,intent(in) :: n 
     double precision,dimension(n),intent(in) :: a,b,c,d 
     double precision,dimension(n),intent(out) :: x 
     double precision,dimension(n) :: cp,dp 
     double precision :: m 
     integer i 

! initialize c-prime and d-prime 
     cp(1) = c(1)/b(1) 
     dp(1) = d(1)/b(1) 
! solve for vectors c-prime and d-prime 
     do i = 2,n 
      m = b(i)-cp(i-1)*a(i) 
      cp(i) = c(i)/m 
      dp(i) = (d(i)-dp(i-1)*a(i))/m 
     enddo 
! initialize x 
     x(n) = dp(n) 
! solve for x from the vectors c-prime and d-prime 
     do i = n-1, 1, -1 
      x(i) = dp(i)-cp(i)*x(i+1) 
     end do 

    end subroutine tridiag 

!Head in terms of moisture 
Function HfTH(humedad) 
    double precision HfTH 
    double precision humedad 
    if (humedad.lt.theta_sat) then 
       HfTH=-1/alpha*(((humedad-theta_res)/(theta_sat-theta_res))**(-1/m)-1)**(1/n) !cm 
    else 
       HfTH=0 
    end if 
    Return 
end function HfTH 

end program EcuacionRichards 

Answer 1

난 당신의 코드에 문제가 임의의 수를 볼 수 있지만 내주의 지속 시간은 아직, 예를 들어, theta_sat를 들어,

당신은 double precision 될 변수의 무리를 선언 단지 가장 악명 그래서 여기에 제한입니다 당신은 디폴트 유형의 리터럴로 그것들을 초기화한다. 문

double precision :: theta_sat=0.43  !cm/cm 

는 0.43double precision 진짜하지 않습니다. 정확히 말하자면, 대부분의 컴파일러에 있을지도 모릅니다. 그리고 컴파일이 기본 실제 변수를 double precision으로 설정하지 않을 때마다 컴파일되지 않습니다. 0.43은 4 바이트 실수이고, theta_sat은 8 바이트 실수이며 컴파일러를 사용하여 theta_sat을 0.43에 가장 가까운 8 바이트 값으로 설정할 수는 없습니다.

현대 Fortran double precision에서는 이전 버전과의 호환성을 유지하지만 더 이상 사용되지 않으므로 kind 유형의 변수 종류를 지정하는 것이 좋습니다. 그래서 어떻게해야하는지에 대한 제안이 가득합니다. 나의 마음에 드는 다음과 같이 고유 모듈 iso_fortran_env에 정의 된 상수를 사용하는 것입니다

use, intrinsic :: iso_fortran_env 

다음과 같이 변수를 선언 :

real(real64) :: theta_sat=0.43_real64  !cm/cm 

참고 값 종류 사양 _real64의 추기.

귀하의 알고리즘이 귀하의 부분에 대한이 실수가 내가 모르는 결과에 실질적으로 영향을 미치기에 충분히 민감한 지 여부.

마지막으로, 프로그램이 올바르지 않다고 말하면 올바르지 않은 방식으로 침묵합니다.