Skip to content

Commit

Permalink
mode debug code
Browse files Browse the repository at this point in the history
  • Loading branch information
@oksanaguba
oksanaguba committed Aug 22, 2019
1 parent 1e85d41 commit 50ffcda
Showing 1 changed file with 86 additions and 15 deletions.
101 changes: 86 additions & 15 deletions components/homme/src/theta-l/prim_advance_mod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -266,6 +266,8 @@ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, co
max_itercnt_perstep = 0
max_itererr_perstep = 0.0

#if 1

a1 = 1/4d0
a2 = 1/6d0
a3 = 3/8d0
Expand Down Expand Up @@ -293,6 +295,7 @@ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, co
itertol=1e-12
call compute_stage_value_dirk(np1,qn0,dhat1*dt,elem,hvcoord,hybrid,&
deriv,nets,nete,maxiter,itertol,tl)

max_itercnt_perstep = max(maxiter,max_itercnt_perstep)
max_itererr_perstep = max(itertol,max_itererr_perstep)

Expand Down Expand Up @@ -324,15 +327,67 @@ subroutine prim_advance_exp(elem, deriv, hvcoord, hybrid,dt, tl, nets, nete, co
deriv,nets,nete,.false.,0d0,1d0,ahat4/a4,1d0)
maxiter=mi
itertol=1e-12

!print *, 'BEFORE LAST DIRK_________________________________________________________'

call compute_stage_value_dirk(np1,qn0,dhat4*dt,elem,hvcoord,hybrid,&
deriv,nets,nete,maxiter,itertol,tl)
max_itercnt_perstep = max(maxiter,max_itercnt_perstep)
max_itererr_perstep = max(itertol,max_itererr_perstep)

call compute_andor_apply_rhs(np1,n0,np1,qn0,a5*dt,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,eta_ave_w,1d0,ahat5/a5,1d0)
#if 1
call compute_andor_apply_rhs(np1,n0,np1,qn0,a5*dt,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,eta_ave_w,1d0,ahat5/a5,1d0)
#else
call compute_andor_apply_rhs(np1,n0,np1,qn0,a5*dt,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,eta_ave_w,1d0,0d0,1d0)
call compute_stage_value_dirk(np1,qn0,dt,elem,hvcoord,hybrid,&
deriv,nets,nete,maxiter,itertol)
#endif

avg_itercnt = ((nstep)*avg_itercnt + max_itercnt_perstep)/(nstep+1)

#endif



#if 0
! Ullrich 3nd order 5 stage: CFL=sqrt( 4^2 -1) = 3.87
! u1 = u0 + dt/5 RHS(u0) (save u1 in timelevel nm1)
call compute_andor_apply_rhs(nm1,n0,n0,qn0,dt/5,elem,hvcoord,hybrid,&
deriv,nets,nete,compute_diagnostics,eta_ave_w/4,1.d0,1.d0,1.d0)
! u2 = u0 + dt/5 RHS(u1)
call compute_andor_apply_rhs(np1,n0,nm1,qn0,dt/5,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,0d0,1.d0,1.d0,1.d0)
! u3 = u0 + dt/3 RHS(u2)
call compute_andor_apply_rhs(np1,n0,np1,qn0,dt/3,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,0d0,1.d0,1.d0,1.d0)
! u4 = u0 + 2dt/3 RHS(u3)
call compute_andor_apply_rhs(np1,n0,np1,qn0,2*dt/3,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,0d0,1.d0,1.d0,1.d0)
! compute (5*u1/4 - u0/4) in timelevel nm1:
do ie=nets,nete
elem(ie)%state%v(:,:,:,:,nm1)= (5*elem(ie)%state%v(:,:,:,:,nm1) &
- elem(ie)%state%v(:,:,:,:,n0) ) /4
elem(ie)%state%vtheta_dp(:,:,:,nm1)=(5*elem(ie)%state%vtheta_dp(:,:,:,nm1) &
- elem(ie)%state%vtheta_dp(:,:,:,n0) )/4
elem(ie)%state%dp3d(:,:,:,nm1)= (5*elem(ie)%state%dp3d(:,:,:,nm1) &
- elem(ie)%state%dp3d(:,:,:,n0) )/4
elem(ie)%state%w_i(:,:,1:nlevp,nm1)=(5*elem(ie)%state%w_i(:,:,1:nlevp,nm1) &
- elem(ie)%state%w_i(:,:,1:nlevp,n0) )/4
elem(ie)%state%phinh_i(:,:,1:nlev,nm1)=(5*elem(ie)%state%phinh_i(:,:,1:nlev,nm1) &
- elem(ie)%state%phinh_i(:,:,1:nlev,n0) )/4
enddo
! u5 = (5*u1/4 - u0/4) + 3dt/4 RHS(u4)
call compute_andor_apply_rhs(np1,nm1,np1,qn0,3*dt/4,elem,hvcoord,hybrid,&
deriv,nets,nete,.false.,3*eta_ave_w/4,1.d0,1.d0,1.d0)


call compute_stage_value_dirk(np1,qn0,dt,elem,hvcoord,hybrid,&
deriv,nets,nete,maxiter,itertol,tl)
#endif


!================================================================================
elseif (tstep_type == 8) then ! IMKG253, might be more efficient than IMKG254, might be a teeny bit bad at coarse resolution

Expand Down Expand Up @@ -2288,6 +2343,11 @@ end subroutine compute_andor_apply_rhs
!===========================================================================================================
!===========================================================================================================

!#define ORIGD
#define DEBUGD


#ifdef ORIGD
!ORIGINAL
subroutine compute_stage_value_dirk(np1,qn0,dt2,elem,hvcoord,hybrid,&
deriv,nets,nete,maxiter,itertol,tl)
Expand Down Expand Up @@ -2368,7 +2428,7 @@ subroutine compute_stage_value_dirk(np1,qn0,dt2,elem,hvcoord,hybrid,&
#endif

!guess 2
#if 0
#if 1
! hydrostatic pressure
pi_i(:,:,1)=hvcoord%hyai(1)*hvcoord%ps0
do k=1,nlev
Expand Down Expand Up @@ -2678,6 +2738,7 @@ subroutine compute_stage_value_dirk(np1,qn0,dt2,elem,hvcoord,hybrid,&

end subroutine compute_stage_value_dirk
!ORIGINAL
#endif



Expand All @@ -2688,10 +2749,9 @@ end subroutine compute_stage_value_dirk



#ifdef DEBUGD



subroutine compute_stage_value_dirk_(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,nete,maxiter,itertol,tl)
subroutine compute_stage_value_dirk(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,nete,maxiter,itertol,tl)

use physical_constants, only : p0, kappa, Rgas

Expand Down Expand Up @@ -2734,7 +2794,7 @@ subroutine compute_stage_value_dirk_(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,

integer :: i,j,k,l,ie,itercount,info(np,np),itercountmax
integer :: nsafe
integer :: location(2)
integer :: location(2), loc(3)

real (kind=real_kind) :: exnerpi(np,np,nlev), dpds(np,np,nlev), pi_i(np,np,nlevp), pi(np,np,nlev), aa

Expand Down Expand Up @@ -2899,7 +2959,7 @@ subroutine compute_stage_value_dirk_(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,
numi0 = itercount-1
#endif

#if 0
#if 1
!!!! GUESS 1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!setup each iteration this way:
Expand Down Expand Up @@ -3041,7 +3101,7 @@ subroutine compute_stage_value_dirk_(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,
#endif


#if 0
#if 1
!!!! GUESS 2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!setup each iteration this way:
Expand Down Expand Up @@ -3186,20 +3246,31 @@ subroutine compute_stage_value_dirk_(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,
!!!! COMPARE GUESSES
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#if 0
#if 1
aa=maxval(abs(guess0 - guess1))
if(aa > 1e-7) then
if(aa > 1e-6) then
print *, 'max|g0-g1| = ', aa, ' at my GID ', elem(ie)%GlobalID
loc = maxloc(abs(guess0 - guess1))
print *, 'VALUES in question are ', guess0(loc(1),loc(2),loc(3)),guess1(loc(1),loc(2),loc(3)), &
guess0(loc(1),loc(2),loc(3))-guess1(loc(1),loc(2),loc(3))
endif

aa=maxval(abs(guess0 - guess2))
if(aa > 1e-7) then
if(aa > 1e-6) then
print *, 'max|g0-g2| = ', aa, ' at my GID ', elem(ie)%GlobalID
loc = maxloc(abs(guess0 - guess2))
print *, 'VALUES in question are ',&
guess0(loc(1),loc(2),loc(3)),guess2(loc(1),loc(2),loc(3)), &
guess0(loc(1),loc(2),loc(3))-guess2(loc(1),loc(2),loc(3))
endif

aa=maxval(abs(guess1 - guess2))
if(aa > 1e-7) then
if(aa > 1e-6) then
print *, 'max|g1-g2| = ', aa, ' at my GID ', elem(ie)%GlobalID
loc = maxloc(abs(guess1 - guess2))
print *, 'VALUES in question are ',&
guess1(loc(1),loc(2),loc(3)),guess2(loc(1),loc(2),loc(3)), &
guess1(loc(1),loc(2),loc(3))-guess2(loc(1),loc(2),loc(3))
endif

!now see if we can get pt where mu is away from 1
Expand Down Expand Up @@ -3255,11 +3326,11 @@ subroutine compute_stage_value_dirk_(np1,qn0,dt2,elem,hvcoord,hybrid,deriv,nets,
itertol=itererrmax
call t_stopf('compute_stage_value_dirk')

end subroutine compute_stage_value_dirk_

end subroutine compute_stage_value_dirk



#endif



Expand Down

0 comments on commit 50ffcda

Please sign in to comment.