[Commits] [svn:einsteintoolkit] NullNews/branches/tapir/src/ (Rev. 24)

bela at caltech.edu bela at caltech.edu
Tue Jan 7 16:37:25 CST 2014


User: szilagyi
Date: 2014/01/07 04:37 PM

Modified:
 /branches/tapir/src/
  NullNews_GetNews.F90

Log:
 work on variable time-step compatibility and on diagnostic code

File Changes:

Directory: /branches/tapir/src/
===============================

File [modified]: NullNews_GetNews.F90
Delta lines: +65 -15
===================================================================
--- branches/tapir/src/NullNews_GetNews.F90	2014-01-07 22:36:50 UTC (rev 23)
+++ branches/tapir/src/NullNews_GetNews.F90	2014-01-07 22:37:24 UTC (rev 24)
@@ -12,6 +12,7 @@
   use NullNews_ScriUtil, only: NullNews_ResetInactiveRe
   use NullInterp 
   use NullDecomp_SpinDecomp, only: SpinDecompFilter
+  use NullDecomp_IO
 
   implicit none
 
@@ -24,20 +25,22 @@
   CCTK_REAL,    dimension (:,:,:), allocatable, save::&
        beta, omega, beta_u, u0
 
-  CCTK_REAL, dimension(3) :: BUF
   logical, save :: FirstTime = .True.
   logical, save :: ThisProc = .false.
   integer, save :: loc_q, loc_p, myproc
 
   CCTK_INT, dimension(1) :: Tarr
-  character(len=500) :: message
 
   CCTK_INT :: ip
+  CCTK_REAL :: dt
 
   DECLARE_CCTK_FUNCTIONS
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_PARAMETERS
 
+  logical :: truncate
+  truncate = (IO_TruncateOutputFiles(cctkGH) .ne. 0) .and. FirstTime
+
   !   call CCTK_INFO("Null News, Calculating News")
   if (FirstTime) then
      FirstTime = .false.
@@ -82,9 +85,16 @@
      endif
      myproc = CCTK_MyProc(cctkGH)     
   endif
+
+  ! assumes 'time' is time on the new level
+  time_of_news = 0.5*(null_time+null_time_p)
+  dt=null_time-null_time_p
+
   if (use_linearized_omega .eq. 0) then
-    call NullNews_integ_omega (cctkGH, null_lsh, cctk_delta_time, tmp_rgfn, tmp_rgfs, Uo, Un, omegao, omegan,betao,betan)
-    call NullNews_integ_comega (cctkGH, null_lsh, cctk_delta_time, tmp_cgfn, tmp_cgfs, omegao, omegan, comegao, comegan)
+    call NullNews_integ_omega (cctkGH, null_lsh, dt, &
+      tmp_rgfn, tmp_rgfs, Uo, Un, omegao, omegan,betao,betan)
+    call NullNews_integ_comega (cctkGH, null_lsh, dt, &
+      tmp_cgfn, tmp_cgfs, omegao, omegan, comegao, comegan)
   else 
     call NullNews_lin_omega (cctkGH, null_lsh, zeta, Jn, omegan, comegan)
   endif
@@ -107,19 +117,56 @@
   ! mid level values, from the characteristic evolution
 
   J = 0.5d0 * (Jn + Jo)
-  J_u = (Jn - Jo) / cctk_delta_time
+  J_u = (Jn - Jo) / dt
   J_l = 0.5d0 * (Jn_l + Jo_l)
-  J_l_u = (Jn_l - Jo_l) / cctk_delta_time
+  J_l_u = (Jn_l - Jo_l) / dt
   beta = 0.5d0 * (betan + betao)
   cB = 0.5d0 * (cBn + cBo)
   U = 0.5d0 * (Un + Uo)
 
   ! For psi4 calculation
-  beta_u = (betan - betao) / cctk_delta_time
-  U_u = (Un - Uo) / cctk_delta_time
+  beta_u = (betan - betao) / dt
+
+  U_u = (Un - Uo) / dt
   U_l = 0.5d0 * (Un_l + Uo_l)
   U_l_l = 0.5d0 * (Un_l_l + Uo_l_l)
 
+
+  if (debug_output == 1) then
+
+      call NullDecomp_WriteCoefsForComplexFunc('J_scri', cctkGH, truncate, &
+            null_lsh, 2, zeta, J(:,:,1), J(:,:,2), time_of_news)
+
+      call NullDecomp_WriteCoefsForComplexFunc('Ju_scri', cctkGH, truncate, &
+            null_lsh, 2, zeta, J_u(:,:,1), J_u(:,:,2), time_of_news)
+
+      call NullDecomp_WriteCoefsForComplexFunc('Jl_scri', cctkGH, truncate, &
+            null_lsh, 2, zeta, J_l(:,:,1), J_l(:,:,2), time_of_news)
+
+      call NullDecomp_WriteCoefsForComplexFunc('Jlu_scri', cctkGH, truncate, &
+            null_lsh, 2, zeta, J_l_u(:,:,1), J_l_u(:,:,2), time_of_news)
+
+     call NullDecomp_WriteCoefsForRealFunc('B_scri', cctkGH, truncate, &
+          null_lsh, zeta, beta(:,:,1), beta(:,:,2), time_of_news)
+
+     call NullDecomp_WriteCoefsForRealFunc('Bu_scri', cctkGH, truncate, &
+          null_lsh, zeta, beta_u(:,:,1), beta_u(:,:,2), time_of_news)
+
+     call NullDecomp_WriteCoefsForComplexFunc('U_scri', cctkGH, truncate, &
+            null_lsh, 1, zeta, U(:,:,1), U(:,:,2), time_of_news)
+
+      call NullDecomp_WriteCoefsForComplexFunc('Uu_scri', cctkGH, truncate, &
+            null_lsh, 1, zeta, U_u(:,:,1), U_u(:,:,2), time_of_news)
+
+      call NullDecomp_WriteCoefsForComplexFunc('Ul_scri', cctkGH, truncate, &
+            null_lsh, 1, zeta, U_l(:,:,1), U_l(:,:,2), time_of_news)
+
+      call NullDecomp_WriteCoefsForComplexFunc('Ull_scri', cctkGH, truncate, &
+            null_lsh, 1, zeta, U_l_l(:,:,1), U_l_l(:,:,2), time_of_news)
+
+  endif
+
+
   call NullInterp_d1 (eth_U(:,:,1), U(:,:,1), 1, 1)
   call NullInterp_d1 (eth_U(:,:,2), U(:,:,2), 1, 1)
 
@@ -154,9 +201,6 @@
   call NullNews_ResetInactive(null_lsh, eth_J)
   call NullNews_ResetInactive(null_lsh, beth_J)
 
-  ! assumes 'time' is time on the new level
-  time_of_news = cctk_time - .5d0 * cctk_delta_time    
-
   ! mid level values, from quantities computed above
 
   sigmaJo = sigmaJn; sigmaKo = sigmaKn; sigmauo = sigmaun
@@ -170,7 +214,7 @@
 
   call NullNews_Psi4 (cctkGH, null_lsh, tmp_cgfn, tmp_cgfs, Jo, Jo_l,betao, cBo, Uo, Uo_l,&
        sigmaJn, sigmaKn, sigmaun,sigmarn, sigmarun, sigmauun,&
-       sigmaJo, sigmaKo, sigmauo,sigmaro, sigmaruo, sigmauuo,Psi4(:,:,1:2),cctk_delta_time)
+       sigmaJo, sigmaKo, sigmauo,sigmaro, sigmaruo, sigmauuo,Psi4(:,:,1:2),dt)
 
   Psi4(:,:,1:2) = Psi4(:,:,1:2)/omega**3
 
@@ -185,7 +229,7 @@
   endif
 
   if (DEBUG_skip_BondiNews.eq.0) &
-       call NullNews_News2BondiNews(cctkGH, null_lsh, cctk_delta_time, qsize,&
+       call NullNews_News2BondiNews(cctkGH, null_lsh, dt, qsize,&
        Uo, Un, Uyo, Uyn, zEvolo, zEvoln,&
        zEvolh, News(:,:,1:2), NewsB(:,:,1:2), uBondio, uBondi, uBondin,&
        redshiftB, omega, J, eth_J, beth_J, eth_U, beth_U,&
@@ -194,8 +238,15 @@
        n_tmp3_cgfn, n_tmp3_cgfs, n_tmp1_rgfn, n_tmp1_rgfs,&
        n_tmp2_rgfn, n_tmp2_rgfs, starting, EV_mask)
 
-  dMdOmega = NewsB(:,:,1:2) * conjg(NewsB(:,:,1:2)) * redshiftB * circle
+  if (debug_output == 1) then
 
+     call NullDecomp_WriteCoefsForRealFunc('uB_scri', cctkGH, truncate, &
+          null_lsh, zeta, uBondi(:,:,1), uBondi(:,:,2), time_of_news)
+
+  end if
+
+  dMdOmega = real(NewsB(:,:,1:2) * conjg(NewsB(:,:,1:2))) * redshiftB * circle
+
   if (mask_NewsB .ne. 0) then
      do ip = 1, 2
         NewsB(:,:,ip) = EQ_mask(:,:) * NewsB(:,:,ip) ! EV_mask
@@ -226,7 +277,6 @@
   omegao = omegan
   comegao = comegan
 
-
 end subroutine NullNews_GetNews
 
 



More information about the Commits mailing list