diff options
| -rw-r--r-- | interface.ccl | 10 | ||||
| -rw-r--r-- | param.ccl | 20 | ||||
| -rw-r--r-- | schedule.ccl | 115 | ||||
| -rw-r--r-- | src/EHFinder_Constants.F90 | 3 | ||||
| -rw-r--r-- | src/EHFinder_SetMask.F90 | 16 | ||||
| -rw-r--r-- | src/EHFinder_SetSym.F90 | 144 | ||||
| -rw-r--r-- | src/EHFinder_Sources2.F90 | 2 | ||||
| -rw-r--r-- | src/EHFinder_Sources3.F90 | 2 | ||||
| -rw-r--r-- | src/EHFinder_mod.F90 | 2 | ||||
| -rw-r--r-- | src/include/centered_second.h | 11 | ||||
| -rw-r--r-- | src/include/physical_part.h | 3 | ||||
| -rw-r--r-- | src/include/upwind_first.h | 3 | ||||
| -rw-r--r-- | src/include/upwind_second.h | 3 | ||||
| -rw-r--r-- | src/make.code.deps | 3 |
14 files changed, 268 insertions, 69 deletions
diff --git a/interface.ccl b/interface.ccl index d543c91..06197d8 100644 --- a/interface.ccl +++ b/interface.ccl @@ -24,13 +24,17 @@ CCTK_REAL dlevel_set type=GF TimeLevels=1 dfx, dfy, dfz, dfsq } "Derivatives of the level set function" -CCTK_REAL ftmp type=GF TimeLevels=1 +CCTK_REAL ftmp_set type=GF TimeLevels=1 +{ + ftmp, sftmp +} "temporary variables used in pde re-parametrization" CCTK_INT eh_mask_all type=GF TimeLevels=1 { eh_mask, tm_mask } "Masks to define active cells" -CCTK_REAL rep_mask type=GF TimeLevels=1 +CCTK_INT rep_mask type=GF TimeLevels=1 -CCTK_INT re_param_control type=SCALAR +CCTK_INT re_param_control_pde type=SCALAR +CCTK_INT re_param_control_approx type=SCALAR @@ -62,14 +62,15 @@ KEYWORD mode "Mode of operation" KEYWORD re_param_method "Integration method in re-parametrization" { "approx" :: "Approximate re-parametrization scheme" - "pde" :: "Re-parametrize by solving an pde" -} "approx" + "pde" :: "Re-parametrize by solving an pde" + "mixed" :: "Use both schemes interchangably" +} "pde" KEYWORD re_param_int_method "Integration method in pde re-parametrization" { "euler" :: "Standard euler scheme" "rk2" :: "Second order Runge-Kutta scheme" -} "rk2" +} "euler" INT re_param_max_iter "maximum number of iteration in the re-parametrization" { @@ -81,12 +82,17 @@ KEYWORD pde_differences "Type of finite diffencing used in pde re-parametrizatio "centered" :: "Use 2nd order centered differences except at the boundaries" "upwind" :: "Use 1st order upwinded differences everywhere" "upwind2" :: "Use 2nd order upwinded differences everywhere" -} "upwind" +} "upwind2" -INT reparametrize_every "Re-parametrize every" +INT reparametrize_every_pde "Re-parametrize every using pde method" { -0: :: "If 0 don't re-parametrize" -} 40 +0: :: "If 0 don't re-parametrize using pde method" +} 100 + +INT reparametrize_every_approx "Re-parametrize every using approx method" +{ +0: :: "If 0 don't re-parametrize using approx method" +} 10 shares: grid diff --git a/schedule.ccl b/schedule.ccl index fa555bd..d2fffb9 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -6,7 +6,7 @@ STORAGE: slevel_set STORAGE: dlevel_set STORAGE: eh_mask_all STORAGE: rep_mask -STORAGE: ftmp +STORAGE: ftmp_set # Check for metric_state @@ -68,24 +68,26 @@ if (CCTK_Equals(mode,"normal")) { } "Re-parametrize the level set function" - if (CCTK_Equals(re_param_method,"pde")) + if (CCTK_Equals(re_param_method,"pde") || CCTK_Equals(re_param_method,"approx") || CCTK_Equals(re_param_method,"mixed")) { - - # Schedule the control routine, that initialises re_param_control if - # it is time for re-parametrization + # Schedule the control routine, that initialises re_param_control_approx + # and re_param_control_pde if it is time for re-parametrization schedule EHFinder_ReParametrizeControl in EHFinder_ReParametrize { LANG: Fortran } "Initializes the re-parametrization control" + } + if (CCTK_Equals(re_param_method,"pde") || CCTK_Equals(re_param_method,"mixed")) + { if (CCTK_Equals(re_param_int_method,"euler")) { # Set up the schedule group for euler re-parametrization using - # re_param_control to control the exit from the while loop. + # re_param_control_pde to control the exit from the while loop. - schedule GROUP Euler_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrize1 WHILE ehfinder::re_param_control + schedule GROUP Euler_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrizeControl WHILE ehfinder::re_param_control_pde { } "Schedule group for Euler re-parametrization" @@ -99,12 +101,95 @@ if (CCTK_Equals(mode,"normal")) # Then synchronize the level set and apply symmetry boundary conditions. - schedule EHFinder_ApplySym in Euler_ReParametrize after EHFinder_ReParametrizeEuler + schedule EHFinder_ApplySymF in Euler_ReParametrize after EHFinder_ReParametrizeEuler { LANG: Fortran SYNC: level_set - } "apply symmetry boundaries" + } "Apply symmetry boundaries and sync" } + + if (CCTK_Equals(re_param_int_method,"rk2")) + { + + # Set up the schedule group for rk2 re-parametrization using + # re_param_control_pde to control the exit from the while loop. + + schedule GROUP RK2_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrize1 WHILE ehfinder::re_param_control_pde + { + } "Schedule group for RK2 re-parametrization" + + # Schedule the routine that does the first RK2 step + + schedule EHFinder_ReParametrizeRK2_1 in RK2_ReParametrize + { + LANG: Fortran + } "Euler scheme" + + # Then synchronize the level set and apply symmetry boundary conditions. + + schedule EHFinder_ApplySymF in RK2_ReParametrize after EHFinder_ReParametrizeRK2_1 + { + LANG: Fortran + SYNC: level_set + } "Apply symmetry boundaries and sync" + + # Schedule the routine that does the second RK2 step + + schedule EHFinder_ReParametrizeRK2_2 in RK2_ReParametrize after EHFinder_ApplySymF + { + LANG: Fortran + } "Euler scheme" + + # Then synchronize the level set and apply symmetry boundary conditions. + + schedule EHFinder_ApplySymF in RK2_ReParametrize after EHFinder_ReParametrizeRK2_2 + { + LANG: Fortran + SYNC: level_set + } "Apply symmetry boundaries and sync" + } + } + + if (CCTK_Equals(re_param_method,"approx") || CCTK_Equals(re_param_method,"mixed")) + { + + # Schedule a routine to update the points closest to the zero-level surface. + + schedule EHFinder_ReParametrize5 in EHFinder_ReParametrize after EHFinder_ReParametrizeControl + { + LANG: Fortran + } "Update the points closest to the surface" + + schedule EHFinder_ApplySymFRep in EHFinder_ReParametrize after EHFinder_ReParametrize5 + { + LANG: Fortran + SYNC: level_set + SYNC: rep_mask + } "Apply symmetry boundaries and sync" + + # Set up the schedule group for approx re-parametrization using + # re_param_control_approx to control the exit from the while loop. + + schedule GROUP Approx_ReParametrize in EHFinder_ReParametrize AFTER EHFinder_ReParametrizeControl WHILE ehfinder::re_param_control_approx + { + } "Schedule group for Approximate re-parametrization" + + schedule EHFinder_ReParametrize6 in Approx_ReParametrize + { + LANG: Fortran + } "Update cells with neighbours that have already been updated" + + schedule EHFinder_ReParametrize7 in Approx_ReParametrize after ReParametrize6 + { + LANG: Fortran + } "Update rep_mask" + + schedule EHFinder_ApplySymFRep in Approx_ReParametrize after EHFinder_ReParametrize7 + { + LANG: Fortran + SYNC: level_set + SYNC: rep_mask + } "Apply symmetry boundaries and sync" } @@ -125,20 +210,26 @@ if (CCTK_Equals(mode,"normal")) # Then apply the symmetry boundary conditions. - schedule EHFinder_ApplySym in EHFinder_SetMask after EHFinder_SetMask1 + schedule EHFinder_ApplySymAll in EHFinder_SetMask after EHFinder_SetMask1 { LANG: Fortran - } "apply symmetry boundaries" + } "Apply symmetry boundaries" # Finally locate the mask boundary and add values to distinguish different # directions. - schedule EHFinder_SetMask2 in EHFinder_SetMask after EHFinder_ApplySym + schedule EHFinder_SetMask2 in EHFinder_SetMask after EHFinder_ApplySymAll { LANG: Fortran SYNC: eh_mask_all } "Finish modifying the mask" + schedule EHFinder_ApplySymMask in EHFinder_SetMask after EHFinder_SetMask2 + { + LANG: Fortran + } "Apply symmetry boundaries" + + } diff --git a/src/EHFinder_Constants.F90 b/src/EHFinder_Constants.F90 index 3da3c8e..fa86316 100644 --- a/src/EHFinder_Constants.F90 +++ b/src/EHFinder_Constants.F90 @@ -1,3 +1,6 @@ +! Various useful constants. +! $Header$ + #include "cctk.h" module EHFinder_Constants diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90 index 995b3a5..6cabfa0 100644 --- a/src/EHFinder_SetMask.F90 +++ b/src/EHFinder_SetMask.F90 @@ -1,4 +1,4 @@ -! Initialisation of the level set function +! Mask modification functions. ! $Header$ #include "cctk.h" @@ -79,8 +79,11 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) active = .false. if ( mask_first ) active = .true. - if ( reparametrize_every .gt. 0 ) then - if ( mod(cctk_iteration,reparametrize_every) .eq. 0 ) active = .true. + if ( reparametrize_every_approx .gt. 0 ) then + if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) active = .true. + end if + if ( reparametrize_every_pde .gt. 0 ) then + if ( mod(cctk_iteration,reparametrize_every_pde) .eq. 0 ) active = .true. end if if ( active ) then @@ -221,8 +224,11 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) active = .false. if ( mask_first ) active = .true. - if ( reparametrize_every .gt. 0 ) then - if ( mod(cctk_iteration,reparametrize_every) .eq. 0 ) active = .true. + if ( reparametrize_every_approx .gt. 0 ) then + if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) active = .true. + end if + if ( reparametrize_every_pde .gt. 0 ) then + if ( mod(cctk_iteration,reparametrize_every_pde) .eq. 0 ) active = .true. end if if ( active ) then diff --git a/src/EHFinder_SetSym.F90 b/src/EHFinder_SetSym.F90 index a8c7660..7559834 100644 --- a/src/EHFinder_SetSym.F90 +++ b/src/EHFinder_SetSym.F90 @@ -23,16 +23,86 @@ subroutine EHFinder_SetSym(CCTK_ARGUMENTS) call CCTK_WARN(1,'Failed to register symmetry for the level set function') end if - call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::rep_mask') - if ( ierr .gt. 0 ) then - call CCTK_WARN(1,'Failed to register symmetry for the level reparametrization mask') - end if +! call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::rep_mask') +! if ( ierr .gt. 0 ) then +! call CCTK_WARN(1,'Failed to register symmetry for the level reparametrization mask') +! end if return end subroutine EHFinder_SetSym -subroutine EHFinder_ApplySym(CCTK_ARGUMENTS) +subroutine EHFinder_ApplySymAll(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + call CartSymVN(ierr,cctkGH,'ehfinder::f') + + # include "include/physical_part.h" + + if ( symx ) then + eh_mask(ngx:1:-1,:,:) = eh_mask(ngx+1:ngx+ngx,:,:) + rep_mask(ngx:1:-1,:,:) = rep_mask(ngx+1:ngx+ngx,:,:) + end if + if ( symy ) then + eh_mask(:,ngy:1:-1,:) = eh_mask(:,ngy+1:ngy+ngy,:) + rep_mask(:,ngy:1:-1,:) = rep_mask(:,ngy+1:ngy+ngy,:) + end if + if ( symz ) then + eh_mask(:,:,ngz:1:-1) = eh_mask(:,:,ngz+1:ngz+ngz) + rep_mask(:,:,ngz:1:-1) = rep_mask(:,:,ngz+1:ngz+ngz) + end if + return +end subroutine EHFinder_ApplySymAll + + +subroutine EHFinder_ApplySymF(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + call CartSymVN(ierr,cctkGH,'ehfinder::f') + return +end subroutine EHFinder_ApplySymF + + +subroutine EHFinder_ApplySymRep(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS + + # include "include/physical_part.h" + + if ( symx ) then + rep_mask(ngx:1:-1,:,:) = rep_mask(ngx+1:ngx+ngx,:,:) + end if + if ( symy ) then + rep_mask(:,ngy:1:-1,:) = rep_mask(:,ngy+1:ngy+ngy,:) + end if + if ( symz ) then + rep_mask(:,:,ngz:1:-1) = rep_mask(:,:,ngz+1:ngz+ngz) + end if + return +end subroutine EHFinder_ApplySymRep + + +subroutine EHFinder_ApplySymFRep(CCTK_ARGUMENTS) use EHFinder_mod @@ -42,37 +112,43 @@ subroutine EHFinder_ApplySym(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS -! print*,'before3' -! print*,rep_mask(24,21,1:4) -! print*,f(24,21,1:4) -! print* + # include "include/physical_part.h" + call CartSymVN(ierr,cctkGH,'ehfinder::f') -! print*,'after3' -! print*,rep_mask(24,21,1:4) -! print*,f(24,21,1:4) -! print* -! -! print*,'before4' -! print*,rep_mask(24,21,1:4) -! print*,f(24,21,1:4) -! print* - call CartSymVN(ierr,cctkGH,'ehfinder::rep_mask') -! -! print*,'after4' -! print*,rep_mask(24,21,1:4) -! print*,f(24,21,1:4) -! print* + + if ( symx ) then + rep_mask(ngx:1:-1,:,:) = rep_mask(ngx+1:ngx+ngx,:,:) + end if + if ( symy ) then + rep_mask(:,ngy:1:-1,:) = rep_mask(:,ngy+1:ngy+ngy,:) + end if + if ( symz ) then + rep_mask(:,:,ngz:1:-1) = rep_mask(:,:,ngz+1:ngz+ngz) + end if + return +end subroutine EHFinder_ApplySymFRep + + +subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS # include "include/physical_part.h" - if ( symx ) then - eh_mask(ngx:1:-1,:,:) = eh_mask(ngx+1:ngx+ngx,:,:) - end if - if ( symy ) then - eh_mask(:,ngy:1:-1,:) = eh_mask(:,ngy+1:ngy+ngy,:) - end if - if ( symz ) then - eh_mask(:,:,ngz:1:-1) = eh_mask(:,:,ngz+1:ngz+ngz) - end if + if ( symx ) then + eh_mask(ngx:1:-1,:,:) = eh_mask(ngx+1:ngx+ngx,:,:) + end if + if ( symy ) then + eh_mask(:,ngy:1:-1,:) = eh_mask(:,ngy+1:ngy+ngy,:) + end if + if ( symz ) then + eh_mask(:,:,ngz:1:-1) = eh_mask(:,:,ngz+1:ngz+ngz) + end if return -end subroutine EHFinder_ApplySym +end subroutine EHFinder_ApplySymMask diff --git a/src/EHFinder_Sources2.F90 b/src/EHFinder_Sources2.F90 index 8a51632..0707317 100644 --- a/src/EHFinder_Sources2.F90 +++ b/src/EHFinder_Sources2.F90 @@ -1,4 +1,4 @@ -! Re-Parametrization of the level set function +! Source function for testing of pde re-parametrization with MoL ! $Header$ #include "cctk.h" diff --git a/src/EHFinder_Sources3.F90 b/src/EHFinder_Sources3.F90 index bb278fc..a0239e5 100644 --- a/src/EHFinder_Sources3.F90 +++ b/src/EHFinder_Sources3.F90 @@ -1,4 +1,4 @@ -! Re-Parametrization of the level set function +! Another (and better) source function for the re-parametrization with MoL. ! $Header$ #include "cctk.h" diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90 index cccc3a3..fbf0e0e 100644 --- a/src/EHFinder_mod.F90 +++ b/src/EHFinder_mod.F90 @@ -1,3 +1,5 @@ +! Module with various common variable. +! $Header$ #include "cctk.h" module EHFinder_mod diff --git a/src/include/centered_second.h b/src/include/centered_second.h index 659e23a..e0d552e 100644 --- a/src/include/centered_second.h +++ b/src/include/centered_second.h @@ -1,10 +1,15 @@ +! Calculation of centered differences. +! $Header$ + +# include "physical_part.h" + idx = half / cctk_delta_space(1) idy = half / cctk_delta_space(2) idz = half / cctk_delta_space(3) -do k = 1, nz - do j = 1, ny - do i = 1, nx +do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr if ( eh_mask(i,j,k) .ge. 0 ) then if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. & ( .not. btest(eh_mask(i,j,k),1) ) ) then diff --git a/src/include/physical_part.h b/src/include/physical_part.h index c14439a..6b5cb21 100644 --- a/src/include/physical_part.h +++ b/src/include/physical_part.h @@ -1,4 +1,5 @@ -! Here I figure out the extension of the physical part of a chunk of a grid. +! Figure out the extent of the physical part of a chunk of a grid. +! $Header$ ! nx, ny and nz contains the full size of the chunk. diff --git a/src/include/upwind_first.h b/src/include/upwind_first.h index 976ccc3..6c99d29 100644 --- a/src/include/upwind_first.h +++ b/src/include/upwind_first.h @@ -1,3 +1,6 @@ +! Calculate first order upwinded differences. +! $Header$ + idx = one / cctk_delta_space(1) idy = one / cctk_delta_space(2) idz = one / cctk_delta_space(3) diff --git a/src/include/upwind_second.h b/src/include/upwind_second.h index 3c33a9a..8ab9ba8 100644 --- a/src/include/upwind_second.h +++ b/src/include/upwind_second.h @@ -1,3 +1,6 @@ +! Calculate second order upwinded differences. +! $Header$ + # include "physical_part.h" idx = half / cctk_delta_space(1) diff --git a/src/make.code.deps b/src/make.code.deps index a7db76a..c27b3a0 100644 --- a/src/make.code.deps +++ b/src/make.code.deps @@ -1,6 +1,5 @@ - -# -*-Makefile-*- # Dependencies file for thorn EHFinder. +# $Header$ EHFinder_mod.F90.o: EHFinder_Constants.F90.o EHFinder_Init.F90.o: EHFinder_mod.F90.o |