Changeset 878 for codes/icosagcm/devel/src

Timestamp:

05/29/19 20:33:00 (5 years ago)

Author:

jisesh

Message:

devel: moved DYSL into compute_caldyn_solver.F90,compute_theta.F90 and compute_NH_geopot.F90

Location:

codes/icosagcm/devel/src

Files:

• dynamics/compute_NH_geopot.F90 (modified) (1 diff)
• dynamics/compute_caldyn_solver.F90 (modified) (1 diff)
• dynamics/compute_theta.F90 (modified) (1 diff)
• kernels_unst/theta.k90 (modified) (1 diff)

codes/icosagcm/devel/src/dynamics/compute_NH_geopot.F90

@@ -6 +6 @@
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
 
-  PUBLIC :: compute_NH_geopot
+#include "../unstructured/unstructured.h90"
+
+  PUBLIC :: compute_NH_geopot,compute_NH_geopot_unst
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(compute_NH_geopot)
+  tau2_g=tau*tau/g
+  g2=g*g
+  gm2 = 1./g2
+  vreff = Treff*cpp/preff*kappa
+  gamma = 1./(1.-kappa)
+
+  BARRIER
+
+  ! compute Phi_star
+  SEQUENCE_C1
+    BODY('1,llm+1')
+      Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau)
+    END_BLOCK
+  END_BLOCK
+
+  ! Newton-Raphson iteration : Phi_il contains current guess value
+  DO iter=1,2 ! 2 iterations should be enough
+    ! Compute pressure, A_ik
+    SELECT CASE(caldyn_thermo)
+    CASE(thermo_theta)
+      {% call() compute_p_and_Aik() %}
+        X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij
+        p_ik(CELL) = preff*(X_ij**gamma)
+        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
+      {% endcall %}
+    CASE(thermo_entropy)
+      {% call() compute_p_and_Aik() %}
+        X_ij = log(vreff*rho_ij) + theta(CELL)/cpp
+        p_ik(CELL) = preff*exp(X_ij*gamma)
+        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
+      {% endcall %}
+    CASE DEFAULT
+        PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
+        STOP
+    END SELECT
+
+    ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
+    ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
+
+    SEQUENCE_C1
+      ! Compute residual R_il and B_il
+      PROLOGUE(1)
+        ! bottom interface l=1
+        ml_g2 = gm2*m_il(CELL)
+        B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot
+        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
+                  + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) )
+      END_BLOCK
+      BODY('2,llm')
+        ! inner interfaces
+        ml_g2 = gm2*m_il(CELL)
+        B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2
+        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
+                  + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL)))
+        ! consistency check : if Wil=0 and initial state is in hydrostatic balance
+        ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2
+        ! and residual = tau^2*(ml+(1/g)dl_pi)=0
+      END_BLOCK
+      EPILOGUE('llm+1')
+        ! top interface l=llm+1
+        ml_g2 = gm2*m_il(CELL)
+        B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2
+        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
+                  + tau2_g*( ptop-p_ik(DOWN(CELL)) )
+      END_BLOCK
+      !
+      ! Forward sweep :
+      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
+      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
+      PROLOGUE(1)
+        X_ij = 1./B_il(CELL)
+        C_ik(CELL) = -A_ik(CELL) * X_ij
+        D_il(CELL) =  R_il(CELL) * X_ij
+      END_BLOCK
+      BODY('2,llm')
+        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
+        C_ik(CELL) = -A_ik(CELL) * X_ij
+        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
+      END_BLOCK
+      EPILOGUE('llm+1')
+        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
+        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
+        ! Back substitution :
+        ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
+        ! + Newton-Raphson update
+        ! top interface l=llm+1
+        x_il(CELL) = D_il(CELL)
+        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
+      END_BLOCK
+      BODY('llm,1,-1')
+        ! Back substitution at lower interfaces
+        x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL))
+        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
+      END_BLOCK
+    END_BLOCK
+
+    IF(debug_hevi_solver) THEN
+       PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
+       PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
+       DO l=1,llm+1
+          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
+       END DO
+       DO l=1,llm+1
+          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))
+       END DO
+    END IF
+
+  END DO ! Newton-Raphson
+
+  BARRIER
+
+  debug_hevi_solver=.FALSE.
+END_BLOCK
+
+#endif END_DYSL
+
+SUBROUTINE compute_NH_geopot_unst(tau, m_ik, m_il, theta, W_il, Phi_il)
+  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+  USE disvert_mod, only : g,Treff,cpp,preff,kappa,caldyn_thermo,thermo_theta, &
+    thermo_entropy
+  USE disvert_mod, ONLY : ptop
+  USE data_unstructured_mod, ONLY : primal_num,edge_num,dual_num,id_NH_geopot,debug_hevi_solver_, &
+    PHI_BOT,pbot,rho_bot,enter_trace, exit_trace
+  FIELD_MASS   :: m_ik, theta   ! IN*2
+  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL
+  NUM :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff
+  INTEGER :: iter
+  LOGICAL :: debug_hevi_solver
+  DECLARE_INDICES
+  NUM :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2
+#define COLUMN 0
+#if COLUMN 
+  NUM1(llm)  :: pk, Ak, Ck
+  NUM1(llm+1):: Rl, Bl, Dl, xl
+#define p_ik(l,ij) pk(l)
+#define A_ik(l,ij) Ak(l)
+#define C_ik(l,ij) Ck(l)
+#define R_il(l,ij) Rl(l)
+#define B_il(l,ij) Bl(l)
+#define D_il(l,ij) Dl(l)
+#define x_il(l,ij) xl(l)
+#else
+  FIELD_MASS :: p_ik, A_ik, C_ik
+  FIELD_GEOPOT :: R_il, B_il, D_il, x_il
+#endif
+
+  debug_hevi_solver=.FALSE.
+!$OMP MASTER
+  debug_hevi_solver = debug_hevi_solver_
+!$OMP END MASTER
+
+#define PHI_BOT(ij) Phi_bot
+  START_TRACE(id_NH_geopot, 7,0,0)
+#include "../kernels_unst/compute_NH_geopot.k90"
+  STOP_TRACE
+
+!$OMP MASTER
+  debug_hevi_solver_ = debug_hevi_solver
+!$OMP END MASTER
+#undef PHI_BOT
+
+#if COLUMN
+#undef p_ik
+#undef A_ik
+#undef C_ik
+#undef R_il
+#undef B_il
+#undef D_il
+#undef x_il
+#endif
+#undef COLUMN
+END SUBROUTINE compute_NH_geopot_unst
 
   SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)

codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90

@@ -4 +4 @@
   PRIVATE
 
+#include "../unstructured/unstructured.h90"
+
   PUBLIC :: compute_caldyn_solver
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(caldyn_solver)
+  !   
+  ! Compute pressure (pres) and Exner function (pk)
+  ! kappa = R/Cp
+  ! 1-kappa = Cv/Cp
+  ! Cp/Cv = 1/(1-kappa)
+  gamma    = 1./(1.-kappa)
+  vreff    = Rd*Treff/preff ! reference specific volume
+  Cvd      = 1./(cpp-Rd)
+  Rd_preff = kappa*cpp/preff
+  FORALL_CELLS_EXT()
+    SELECT CASE(caldyn_thermo)
+    CASE(thermo_theta)
+      ON_PRIMAL
+        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
+        rho_ij = rho_ij*g*rhodz(CELL)
+        X_ij = Rd_preff*theta(CELL,1)*rho_ij
+        ! kappa.theta.rho = p/exner
+        ! => X = (p/p0)/(exner/Cp)
+        !      = (p/p0)^(1-kappa)
+        pres(CELL) = preff*(X_ij**gamma)  ! pressure
+        ! Compute Exner function (needed by compute_caldyn_fast)
+        ! other formulae possible if exponentiation is slow
+        pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner
+      END_BLOCK
+    CASE(thermo_entropy)
+      ON_PRIMAL
+        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
+        rho_ij = rho_ij*g*rhodz(CELL)
+        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd )
+        pres(CELL) = rho_ij*Rd*T_ij
+        pk(CELL)   = T_ij
+      END_BLOCK
+    CASE DEFAULT
+      STOP
+    END SELECT
+  END_BLOCK
+
+  ! We need a barrier here because we compute pres above and do a vertical difference below
+  BARRIER
+
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ON_PRIMAL
+      CST_IFTHEN(IS_BOTTOM_LEVEL)
+        ! Lower BC
+        dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL)
+      CST_ELSEIF(IS_TOP_INTERFACE)
+        ! Top BC
+        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL)
+      CST_ELSE
+        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL)
+      CST_ENDIF
+      W(CELL)    = W(CELL)+tau*dW(CELL) ! update W
+      dPhi(CELL) = g*g*W(CELL)/m_il(CELL)
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we update W above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ON_PRIMAL
+      ! compute du = -0.5*g^2.grad(w^2)
+      berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 )
+    END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT()
+    ON_EDGES
+      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
+    END_BLOCK
+  END_BLOCK
+END_BLOCK
+
+#endif END_DYSL
+
+SUBROUTINE compute_caldyn_solver_unst(tau,rhodz,theta, berni,pres,m_il, pk,geopot,W,dPhi,dW,du)
+  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+  USE earth_const
+  USE trace
+  USE grid_param, ONLY : nqdyn
+  USE disvert_mod, ONLY : ptop
+  USE data_unstructured_mod, ONLY : id_solver,primal_num,dual_num,edge_num,left, right,PHI_BOT, &
+    enter_trace, exit_trace
+  USE compute_NH_geopot_mod, ONLY : compute_NH_geopot_unst
+  NUM, INTENT(IN) :: tau 
+  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, BUF*2
+  FIELD_THETA  :: theta                  ! IN
+  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, BUF
+  FIELD_U      :: du                     ! OUT
+  DECLARE_INDICES
+  NUM :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff
+  REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6  ! FIXME
+#define PHI_BOT(ij) Phi_bot
+#include "../kernels_unst/caldyn_mil.k90"
+  IF(tau>0) THEN ! solve implicit problem for geopotential
+    CALL compute_NH_geopot_unst(tau, rhodz, m_il, theta, W, geopot)
+  END IF
+  START_TRACE(id_solver, 7,0,1)
+#include "../kernels_unst/caldyn_solver.k90"
+  STOP_TRACE
+#undef PHI_BOT
+END SUBROUTINE compute_caldyn_solver_unst
 
   SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)

codes/icosagcm/devel/src/dynamics/compute_theta.F90

@@ -4 +4 @@
   PRIVATE
 
+#include "../unstructured/unstructured.h90"
+
   PUBLIC :: compute_theta
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(theta)
+  IF(caldyn_eta==eta_mass) THEN ! Compute mass
+    ! FIXME : here mass_col is computed from rhodz
+    ! so that the DOFs are the same whatever caldyn_eta
+    ! in DYNAMICO mass_col is prognosed rather than rhodz
+    SEQUENCE_C1
+      PROLOGUE(0)
+        mass_col(HIDX(CELL))=0.
+      END_BLOCK
+      BODY('1,llm')
+          mass_col(HIDX(CELL)) = mass_col(HIDX(CELL)) + rhodz(CELL)
+      END_BLOCK
+    END_BLOCK
+    FORALL_CELLS_EXT()
+      ON_PRIMAL
+        ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on
+        ! pressure rather than mass
+        !        m = mass_dak(CELL)+(mass_col(HIDX(CELL))*g+ptop)*mass_dbk(CELL)
+        !        rhodz(CELL) = m/g 
+        rhodz(CELL) = MASS_DAK(CELL) + mass_col(HIDX(CELL))*MASS_DBK(CELL)
+      END_BLOCK
+    END_BLOCK
+  END IF
+  DO iq=1,nqdyn
+    FORALL_CELLS_EXT()
+      ON_PRIMAL
+        theta(CELL,iq) = theta_rhodz(CELL,iq)/rhodz(CELL)
+      END_BLOCK
+    END_BLOCK
+  END DO
+END_BLOCK
+
+#endif END_DYSL
+
+SUBROUTINE compute_theta_unst(mass_col,rhodz,theta_rhodz, theta)
+  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+  USE data_unstructured_mod, ONLY : id_theta,primal_num,dual_num,edge_num, &
+    enter_trace, exit_trace
+  USE grid_param, ONLY : nqdyn
+  USE disvert_mod, ONLY : mass_dak, mass_dbk, caldyn_eta, eta_mass, ptop
+  FIELD_PS :: mass_col
+  FIELD_MASS :: rhodz
+  FIELD_THETA :: theta, theta_rhodz
+  DECLARE_INDICES
+  NUM :: m
+  START_TRACE(id_theta, 3,0,0) ! primal, dual, edge  
+#define MASS_DAK(l,ij) mass_dak(l)
+#define MASS_DBK(l,ij) mass_dbk(l)
+#include "../kernels_unst/theta.k90"
+#undef MASS_DAK
+#undef MASS_DBK
+  STOP_TRACE
+END SUBROUTINE
 
   SUBROUTINE compute_theta(ps,theta_rhodz, rhodz,theta)

codes/icosagcm/devel/src/kernels_unst/theta.k90

@@ -20 +20 @@
             ! m = mass_dak(l,ij)+(mass_col(ij)*g+ptop)*mass_dbk(l,ij)
             ! rhodz(l,ij) = m/g
-            rhodz(l,ij) = mass_dak(l,ij) + mass_col(ij)*mass_dbk(l,ij)
+            rhodz(l,ij) = MASS_DAK(l,ij) + mass_col(ij)*MASS_DBK(l,ij)
          END DO
       END DO