From 3f867f25271a2e65e8d583ad0a366db43f4bea41 Mon Sep 17 00:00:00 2001 From: Loren Date: Tue, 18 Jun 2024 16:31:02 -0700 Subject: [PATCH 1/8] tried pulling in all my changes from 3 years ago --- src/Physics/Checkpointing.F90 | 245 ++++++++++++++++++++++++---------- 1 file changed, 177 insertions(+), 68 deletions(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index ec308eba..024aff18 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -242,19 +242,22 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Real*8, Intent(InOut) :: fields(:,:,:,:), abterms(:,:,:,:) Integer :: n_r_old, l_max_old, grid_type_old Integer :: i, ierr, mp, lb,ub, ab_offset - Integer :: old_pars(7), fcount(3,2), version + Integer :: old_pars(7 + nsubmax), fcount(3,2), version Integer :: last_iter, last_auto, endian_tag, funit Integer*8 :: found_bytes, expected_bytes, n_r_old_big, l_max_old_big Integer :: read_magnetism = 0, read_hydro = 0 Integer, Allocatable :: rinds(:), gpars(:,:), read_var(:) Real*8 :: dt_pars(3),dt,new_dt - Real*8, Allocatable :: old_radius(:), radius_old(:) + Real*8, Allocatable :: old_radius(:), radius_old_loc(:) Real*8, Allocatable :: tempfield1(:,:,:,:), tempfield2(:,:,:,:) Character*120 :: autostring, dstring, iterstring, access_type Character*256 :: grid_file, checkfile Character*1 :: under_slash Character*13 :: szstr Logical :: legacy_format, fexist + Integer :: ncheby_old(1:nsubmax) + Integer :: idom, n_r_loc, n_r_old_loc, irmin, irmax, irmin_old, irmax_old + Logical :: first_time_interpolating = .True. read_hydro = read_pars(1) read_magnetism = read_pars(2) @@ -378,6 +381,16 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif Endif + ! Get the old ncheby + If (ndomains .gt. 1) Then ! get the old ncheby from Checkpoint + ! main_input + input_file = TRIM(my_path)//trim(checkpoint_prefix)//'/main_input' + Call Get_Old_Ncheby(input_file, ncheby_old) + Else ! for ndomains = 1, user need not specify ncheby + ncheby_old(1) = n_r_old + Endif + + If (ierr .eq. 0) Then ! Verify that all checkpoint files exist and have the correct size. n_r_old_big = n_r_old @@ -409,6 +422,7 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) old_pars(2) = grid_type_old old_pars(3) = l_max_old old_pars(7) = ierr + old_pars(8:8+nsubmax-1) = ncheby_old(1:nsubmax) dt_pars(1) = dt dt_pars(2) = new_dt dt_pars(3) = checkpoint_time @@ -416,9 +430,9 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif If (my_row_rank .eq. 0) Then !2-D broadcast pattern - Call MPI_Bcast(old_pars,7, MPI_INTEGER, 0, pfi%ccomm%comm, ierr) + Call MPI_Bcast(old_pars,7+nsubmax, MPI_INTEGER, 0, pfi%ccomm%comm, ierr) Endif - Call MPI_Bcast(old_pars,7, MPI_INTEGER, 0, pfi%rcomm%comm, ierr) + Call MPI_Bcast(old_pars,7+nsubmax, MPI_INTEGER, 0, pfi%rcomm%comm, ierr) n_r_old = old_pars(1) grid_type_old = old_pars(2) @@ -431,6 +445,8 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) legacy_format=.true. Endif + ncheby_old(1:nsubmax) = old_pars(8:8+nsubmax-1) + If (old_pars(7) .ne. 0) Then ! Something is wrong with this checkpoint. If (my_rank .eq. 0) Then @@ -566,78 +582,107 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif - ! NOW, if n_r_old and grid_type_old are the same, we can copy chtkmp%p1b into abterms and - ! fields. Otherwise, we need to interpolate onto the current grid - ! When we change the checkpointing format, should also store AB terms in cheby-space - If ((n_r_old .ne. n_r) .or. (grid_type_old .ne. grid_type) ) Then - ! Interpolate - ! We will assume the user kept the same radial domain bounds. - ! If they have not, this will end badly. - If (my_rank .eq. 0) Then - Call stdout%print(' ') - Call stdout%print('------ Radial grid has changed.') - Call stdout%print('------ Interpolating onto new grid.') - Write(szstr,'(i13)')grid_type_old - Call stdout%print('------ Old grid_type: '//TRIM(szstr)) - Write(szstr,'(i13)')grid_type - Call stdout%print('------ Current grid_type: '//TRIM(szstr)) - Write(szstr,'(i13)')n_r_old - Call stdout%print('------ Old N_R: '//TRIM(szstr)) - Write(szstr,'(i13)')n_r - Call stdout%print('------ Current N_R: '//TRIM(szstr)) - Call stdout%print(' ') + ! Loop over the domains to set the Chebyshev coefficients, possibly + ! interpolating in radius + irmin_old = n_r_old + irmax_old = n_r_old - ncheby_old(1) + 1 + irmin = n_r + irmax = n_r - ncheby(1) + 1 + Do idom = 1, ndomains + ! NOW, if n_r_old and grid_type_old are the same, we can copy chtkmp%p1b into abterms and + ! fields. Otherwise, we need to interpolate onto the current grid + ! When we change the checkpointing format, should also store AB terms in cheby-space + + ! Loop from inner domain to outer domain + n_r_old_loc = ncheby_old(idom) + n_r_loc = ncheby(idom) + + ! Increment (decrement?) the rmin/rmax indices for each domain switch + If (idom .gt. 1) Then + irmin_old = irmin_old - ncheby_old(idom - 1) + irmax_old = irmax_old - ncheby_old(idom) + irmin = irmin - ncheby(idom - 1) + irmax = irmax - ncheby(idom) Endif - If (n_r_old .lt. n_r) Then - - ! The fields are OK - they are already in chebyshev space - fields(:,:,:,1:numfields) = chktmp%p1b(:,:,:,1:numfields) - - ! The AB terms are stored in physical space (in radius). - ! They need to be transformed, coefficients copied, and transformed back.. - ! First, we need to initialize the old chebyshev grid. - Allocate(radius_old(1:n_r_old)) - Call cheby_info%Init(radius_old,rmin,rmax) ! We assume that rmax and rmin do not change - fcount(:,:) = numfields - Call chktmp2%init(field_count = fcount, config = 'p1a') - Call chktmp2%construct('p1a') - chktmp2%p1a(:,:,:,:) = 0.0d0 - ! Allocate tempfield1, tempfield2 - lb = lbound(chktmp%p1b,3) - ub = ubound(chktmp%p1b,3) - Allocate(tempfield1(1:n_r_old,1:2,lb:ub,1)) - Allocate(tempfield2(1:n_r_old,1:2,lb:ub,1)) - - Do i = 1, numfields - tempfield1(:,:,:,:) = 0.0d0 - tempfield2(:,:,:,:) = 0.0d0 - tempfield1(1:n_r_old,:,:,1) = chktmp%p1b(1:n_r_old,:,:,numfields+i) - call cheby_info%tospec4d(tempfield1,tempfield2) - chktmp2%p1a(1:n_r_old,:,:,i) = tempfield2(1:n_r_old,:,:,1) - Enddo - DeAllocate(tempfield1,tempfield2) + ! Interpolate if necessary (in each domain separately) + If ((n_r_old_loc .ne. n_r_loc) .or. (grid_type_old .ne. grid_type) ) Then + ! Interpolate + ! We will assume the user kept the same radial domain bounds. + ! If they have not, this will end badly. + If (my_rank .eq. 0) Then + Call stdout%print(' ') + If (ndomains .gt. 1) Then ! if more than one domain, specify + ! which one we are in + Write(szstr,'(i13)')idom + Call stdout%print('------ In domain: '//TRIM(szstr)) + Endif + Call stdout%print('------ Radial grid has changed.') + Call stdout%print('------ Interpolating onto new grid.') + Write(szstr,'(i13)')grid_type_old + Call stdout%print('------ Old grid_type: '//TRIM(szstr)) + Write(szstr,'(i13)')grid_type + Call stdout%print('------ Current grid_type: '//TRIM(szstr)) + Write(szstr,'(i13)')n_r_old_loc + Call stdout%print('------ Old N_R: '//TRIM(szstr)) + Write(szstr,'(i13)')n_r_loc + Call stdout%print('------ Current N_R: '//TRIM(szstr)) + Call stdout%print(' ') + Endif + If (n_r_old_loc .lt. n_r_loc) Then + ! The fields are OK - they are already in chebyshev space + fields(irmax:irmax+n_r_old_loc-1,:,:,1:numfields) = chktmp%p1b(irmax_old:irmin_old,:,:,1:numfields) + + ! The AB terms are stored in physical space (in radius). + ! They need to be transformed, coefficients copied, and transformed back.. + ! First, we need to initialize the old chebyshev grid. + Allocate(radius_old_loc(1:n_r_old_loc)) + Call cheby_info%Init(radius_old_loc,radius(irmin),radius(irmax)) ! We assume that the domain bounds do not change + fcount(:,:) = numfields + If (first_time_interpolating) Then ! Can only initialize chktmp2 once + Call chktmp2%init(field_count = fcount, config = 'p1a') + first_time_interpolating = .False. + Endif + Call chktmp2%construct('p1a') + chktmp2%p1a(:,:,:,:) = 0.0d0 + ! Allocate tempfield1, tempfield2 + lb = lbound(chktmp%p1b,3) + ub = ubound(chktmp%p1b,3) + Allocate(tempfield1(1:n_r_old_loc,1:2,lb:ub,1)) + Allocate(tempfield2(1:n_r_old_loc,1:2,lb:ub,1)) + + Do i = 1, numfields + tempfield1(:,:,:,:) = 0.0d0 + tempfield2(:,:,:,:) = 0.0d0 + tempfield1(:,:,:,1) = chktmp%p1b(irmax_old:irmin_old,:,:,numfields+i) + Call cheby_info%tospec4d(tempfield1,tempfield2) + chktmp2%p1a(irmax:irmax+n_r_old_loc-1,:,:,i) = tempfield2(:,:,:,1) + Enddo + DeAllocate(tempfield1,tempfield2) + + + Call chktmp2%construct('p1b') + !Normal transform(p1a,p1b) + Call gridcp%From_Spectral(chktmp2%p1a,chktmp2%p1b) + + abterms(irmax:irmin,:,:,1:numfields) = chktmp2%p1b(irmax:irmin,:,:,1:numfields) + Call cheby_info%destroy() + Call chktmp2%deconstruct('p1a') + Call chktmp2%deconstruct('p1b') + Deallocate(radius_old_loc) + Endif - Call chktmp2%construct('p1b') - !Normal transform(p1a,p1b) - Call gridcp%From_Spectral(chktmp2%p1a,chktmp2%p1b) + Else ! n_r_old_loc = n_r_loc + ! Interpolation is complete, now we just copy into the other arrays + fields(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,1:numfields) + abterms(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,numfields+1:numfields*2) - abterms(:,:,:,1:numfields) = chktmp2%p1b(:,:,:,1:numfields) - Call cheby_info%destroy() - Call chktmp2%deconstruct('p1a') - Call chktmp2%deconstruct('p1b') - Deallocate(radius_old) Endif + Enddo - Else - - ! Interpolation is complete, now we just copy into the other arrays - fields(:,:,:,1:numfields) = chktmp%p1b(:,:,:,1:numfields) - abterms(:,:,:,1:numfields) = chktmp%p1b(:,:,:,numfields+1:numfields*2) - - Endif - Call chktmp%deconstruct('p1b') DeAllocate(old_radius) + Call chktmp%deconstruct('p1b') End Subroutine Read_Checkpoint @@ -697,4 +742,68 @@ Subroutine Load_BC_Mask(bvals) bvals(:,:,:,:) = boundary_mask(:,:,:,:) End Subroutine Load_BC_Mask + Subroutine Get_Old_Ncheby(input_file, ncheby_old) + Implicit None + Character*120, Intent(In) :: input_file + Integer, Intent(Out) :: ncheby_old(1:nsubmax) + Integer :: n_r_tmp, n_theta_tmp, nprow_tmp, npcol_tmp + Real*8 :: rmin_tmp, rmax_tmp + Integer :: npout_tmp, precise_bounds_tmp, grid_type_tmp, l_max_tmp, n_l_tmp + Real*8 :: aspect_ratio_tmp, shell_depth_tmp + !Integer, Parameter :: nsubmax = 256 + Integer :: ncheby_tmp(1:nsubmax) + Real*8 :: domain_bounds_tmp(1:nsubmax+1) + Integer :: dealias_by_tmp(1:nsubmax) + Integer :: n_uniform_domains_tmp + Logical :: uniform_bounds_tmp + + ! Don't overwrite the old problemsize namelist! + n_r_tmp = n_r + n_theta_tmp = n_theta + nprow_tmp = nprow + npcol_tmp = npcol + rmin_tmp = rmin + rmax_tmp = rmax + npout_tmp = npout + precise_bounds_tmp = precise_bounds + grid_type_tmp = grid_type + l_max_tmp = l_max + n_l_tmp = n_l + aspect_ratio_tmp = aspect_ratio + shell_depth_tmp = shell_depth + ncheby_tmp(1:nsubmax) = ncheby(1:nsubmax) + domain_bounds_tmp(1:nsubmax+1) = domain_bounds(1:nsubmax+1) + dealias_by_tmp(1:nsubmax) = dealias_by(1:nsubmax) + n_uniform_domains_tmp = n_uniform_domains + uniform_bounds_tmp = uniform_bounds + + ! Read problemsize_namelist of input_file + Open(unit=20, file=input_file, status="old", position="rewind") + Read(unit=20, nml=problemsize_namelist) + Close(20) + + ! "ncheby_old" is now in ncheby + ncheby_old(1:nsubmax) = ncheby(1:nsubmax) + + ! Reset the problemsize_namelist + n_r = n_r_tmp + n_theta = n_theta_tmp + nprow = nprow_tmp + npcol = npcol_tmp + rmin = rmin_tmp + rmax = rmax_tmp + npout = npout_tmp + precise_bounds = precise_bounds_tmp + grid_type = grid_type_tmp + l_max = l_max_tmp + n_l = n_l_tmp + aspect_ratio = aspect_ratio_tmp + shell_depth = shell_depth_tmp + ncheby(1:nsubmax) = ncheby_tmp(1:nsubmax) + domain_bounds(1:nsubmax+1) = domain_bounds_tmp(1:nsubmax+1) + dealias_by(1:nsubmax) = dealias_by_tmp(1:nsubmax) + n_uniform_domains = n_uniform_domains_tmp + uniform_bounds = uniform_bounds_tmp + End Subroutine Get_Old_Ncheby + End Module Checkpointing From c662f59c87a5f1d9ac1aea4a086e9973caf51d5c Mon Sep 17 00:00:00 2001 From: Loren Date: Tue, 18 Jun 2024 16:42:10 -0700 Subject: [PATCH 2/8] forgot to define input_file type --- src/Physics/Checkpointing.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index 024aff18..1bdafb0c 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -251,7 +251,7 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Real*8, Allocatable :: old_radius(:), radius_old_loc(:) Real*8, Allocatable :: tempfield1(:,:,:,:), tempfield2(:,:,:,:) Character*120 :: autostring, dstring, iterstring, access_type - Character*256 :: grid_file, checkfile + Character*256 :: grid_file, checkfile, input_file Character*1 :: under_slash Character*13 :: szstr Logical :: legacy_format, fexist From 5fad63ae57e7186101c2b50df05966a9d129151f Mon Sep 17 00:00:00 2001 From: Loren Date: Wed, 19 Jun 2024 15:26:27 -0700 Subject: [PATCH 3/8] making comments clearer --- src/Physics/Checkpointing.F90 | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index 1bdafb0c..c9f87972 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -582,8 +582,11 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif - ! Loop over the domains to set the Chebyshev coefficients, possibly - ! interpolating in radius + ! Loop over the domains to set the Chebyshev coefficients, + ! possibly interpolating in radius for each subdomain + ! The fields are easy, since they are stored in spectral (Chebyshev) space + ! The AB terms are a bit more complicated, since they are currently stored in physical (radius) space + ! When we change the checkpointing format, should also store AB terms in cheby-space irmin_old = n_r_old irmax_old = n_r_old - ncheby_old(1) + 1 irmin = n_r @@ -591,13 +594,17 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Do idom = 1, ndomains ! NOW, if n_r_old and grid_type_old are the same, we can copy chtkmp%p1b into abterms and ! fields. Otherwise, we need to interpolate onto the current grid - ! When we change the checkpointing format, should also store AB terms in cheby-space + ! And by "interpolate," we mean set the lowest-order Chebyshev coefficients of + ! the new checkpoint using the old checkpoint, and zero out the higher-order coefficients. + ! For abterms, we will transform to spectral space, set the lowest-order coefficients and + ! and transform back ! Loop from inner domain to outer domain n_r_old_loc = ncheby_old(idom) n_r_loc = ncheby(idom) - ! Increment (decrement?) the rmin/rmax indices for each domain switch + ! Decrement the rmin/rmax indices (which go down with increasing radius/idom + ! for each subdomain after the first If (idom .gt. 1) Then irmin_old = irmin_old - ncheby_old(idom - 1) irmax_old = irmax_old - ncheby_old(idom) @@ -632,10 +639,11 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) If (n_r_old_loc .lt. n_r_loc) Then ! The fields are OK - they are already in chebyshev space + ! Just set the lowest-order coefficients from the checkpoint fields(irmax:irmax+n_r_old_loc-1,:,:,1:numfields) = chktmp%p1b(irmax_old:irmin_old,:,:,1:numfields) ! The AB terms are stored in physical space (in radius). - ! They need to be transformed, coefficients copied, and transformed back.. + ! They need to be transformed, lowest-order coefficients copied, and transformed back. ! First, we need to initialize the old chebyshev grid. Allocate(radius_old_loc(1:n_r_old_loc)) Call cheby_info%Init(radius_old_loc,radius(irmin),radius(irmax)) ! We assume that the domain bounds do not change From 64b3b9504884e863e2652045f3c60e2eb1075c47 Mon Sep 17 00:00:00 2001 From: Loren Date: Wed, 19 Jun 2024 16:16:28 -0700 Subject: [PATCH 4/8] for clarity, added many comments to radial interpolation sequence --- src/Physics/Checkpointing.F90 | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index c9f87972..ab1d4af8 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -660,7 +660,11 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Allocate(tempfield1(1:n_r_old_loc,1:2,lb:ub,1)) Allocate(tempfield2(1:n_r_old_loc,1:2,lb:ub,1)) - Do i = 1, numfields + ! For each field, put the physical AB terms from the checkpoint in tempfield1 + ! use "tospec4d" to put the spectral tempfield1 coefficients in tempfield2 + ! copy the tempfield2 spectral coefficients into chktmp2 + ! (thus interpolating, since chktmp2 will have zeros for the higher-order coefficients) + Do i = 1, numfields tempfield1(:,:,:,:) = 0.0d0 tempfield2(:,:,:,:) = 0.0d0 tempfield1(:,:,:,1) = chktmp%p1b(irmax_old:irmin_old,:,:,numfields+i) @@ -669,12 +673,14 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Enddo DeAllocate(tempfield1,tempfield2) - + ! chktmp2 now has the interpolated abterms from the checkpoint in spectral space + ! use "From_Spectral" to transform these fields into physical space and copy into abterms Call chktmp2%construct('p1b') !Normal transform(p1a,p1b) Call gridcp%From_Spectral(chktmp2%p1a,chktmp2%p1b) - abterms(irmax:irmin,:,:,1:numfields) = chktmp2%p1b(irmax:irmin,:,:,1:numfields) + + ! clean up Call cheby_info%destroy() Call chktmp2%deconstruct('p1a') Call chktmp2%deconstruct('p1b') @@ -682,7 +688,7 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif Else ! n_r_old_loc = n_r_loc - ! Interpolation is complete, now we just copy into the other arrays + ! no interpolation is needed, we just copy the checkpoint into fields and abterms fields(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,1:numfields) abterms(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,numfields+1:numfields*2) From afc2cd8163ba705850cae366cf0a2e30ca4392b3 Mon Sep 17 00:00:00 2001 From: Loren Date: Wed, 19 Jun 2024 20:14:07 -0700 Subject: [PATCH 5/8] added error message and exit if user goes down in resolution --- src/Physics/Checkpointing.F90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index ab1d4af8..c88615b7 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -685,8 +685,15 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Call chktmp2%deconstruct('p1a') Call chktmp2%deconstruct('p1b') Deallocate(radius_old_loc) + Else ! Rayleigh doesn't currently support degrading radial resolution--exit now + If (my_rank .eq. 0) Then + Call stdout%print('ERROR: Rayleigh currently does not support degrading radial resolution.') + Call stdout%print('Now exiting') + Call stdout%partial_flush() + Endif + Call pfi%exit() + Stop Endif - Else ! n_r_old_loc = n_r_loc ! no interpolation is needed, we just copy the checkpoint into fields and abterms fields(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,1:numfields) From 209cd78a57525975972561d32c54da8d87a4ecf7 Mon Sep 17 00:00:00 2001 From: Loren Date: Thu, 20 Jun 2024 11:42:57 -0700 Subject: [PATCH 6/8] made get_old_ncheby() much simpler and less dangerous --- src/Physics/Checkpointing.F90 | 90 +++++++++++------------------------ 1 file changed, 29 insertions(+), 61 deletions(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index c88615b7..db06b6c0 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -385,7 +385,7 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) If (ndomains .gt. 1) Then ! get the old ncheby from Checkpoint ! main_input input_file = TRIM(my_path)//trim(checkpoint_prefix)//'/main_input' - Call Get_Old_Ncheby(input_file, ncheby_old) + Call Get_Old_Ncheby(old_radius, n_r_old, ncheby_old) Else ! for ndomains = 1, user need not specify ncheby ncheby_old(1) = n_r_old Endif @@ -763,68 +763,36 @@ Subroutine Load_BC_Mask(bvals) bvals(:,:,:,:) = boundary_mask(:,:,:,:) End Subroutine Load_BC_Mask - Subroutine Get_Old_Ncheby(input_file, ncheby_old) + Subroutine Get_Old_Ncheby(old_radius, n_r_old, ncheby_old) Implicit None - Character*120, Intent(In) :: input_file + ! Get the old ncheby from old_radius, looking for repeated radial points Integer, Intent(Out) :: ncheby_old(1:nsubmax) - Integer :: n_r_tmp, n_theta_tmp, nprow_tmp, npcol_tmp - Real*8 :: rmin_tmp, rmax_tmp - Integer :: npout_tmp, precise_bounds_tmp, grid_type_tmp, l_max_tmp, n_l_tmp - Real*8 :: aspect_ratio_tmp, shell_depth_tmp - !Integer, Parameter :: nsubmax = 256 - Integer :: ncheby_tmp(1:nsubmax) - Real*8 :: domain_bounds_tmp(1:nsubmax+1) - Integer :: dealias_by_tmp(1:nsubmax) - Integer :: n_uniform_domains_tmp - Logical :: uniform_bounds_tmp - - ! Don't overwrite the old problemsize namelist! - n_r_tmp = n_r - n_theta_tmp = n_theta - nprow_tmp = nprow - npcol_tmp = npcol - rmin_tmp = rmin - rmax_tmp = rmax - npout_tmp = npout - precise_bounds_tmp = precise_bounds - grid_type_tmp = grid_type - l_max_tmp = l_max - n_l_tmp = n_l - aspect_ratio_tmp = aspect_ratio - shell_depth_tmp = shell_depth - ncheby_tmp(1:nsubmax) = ncheby(1:nsubmax) - domain_bounds_tmp(1:nsubmax+1) = domain_bounds(1:nsubmax+1) - dealias_by_tmp(1:nsubmax) = dealias_by(1:nsubmax) - n_uniform_domains_tmp = n_uniform_domains - uniform_bounds_tmp = uniform_bounds - - ! Read problemsize_namelist of input_file - Open(unit=20, file=input_file, status="old", position="rewind") - Read(unit=20, nml=problemsize_namelist) - Close(20) - - ! "ncheby_old" is now in ncheby - ncheby_old(1:nsubmax) = ncheby(1:nsubmax) - - ! Reset the problemsize_namelist - n_r = n_r_tmp - n_theta = n_theta_tmp - nprow = nprow_tmp - npcol = npcol_tmp - rmin = rmin_tmp - rmax = rmax_tmp - npout = npout_tmp - precise_bounds = precise_bounds_tmp - grid_type = grid_type_tmp - l_max = l_max_tmp - n_l = n_l_tmp - aspect_ratio = aspect_ratio_tmp - shell_depth = shell_depth_tmp - ncheby(1:nsubmax) = ncheby_tmp(1:nsubmax) - domain_bounds(1:nsubmax+1) = domain_bounds_tmp(1:nsubmax+1) - dealias_by(1:nsubmax) = dealias_by_tmp(1:nsubmax) - n_uniform_domains = n_uniform_domains_tmp - uniform_bounds = uniform_bounds_tmp + Integer, Intent(In) :: n_r_old + Real*8, Intent(In) :: old_radius(:) + Real*8 :: tol = 1.0d-10 ! don't check if things are exactly equal + Real*8 :: r_loc, r_above + Integer :: ir, nr_in_lower_domains, idom + + ! Loop from bottom to top of shell + ! Check if the grid point just above (index - 1) is the same (to within relative tol) + ! If so, the current index marks the top of a subdomain + nr_in_lower_domains = 0 + idom = 1 + Do ir = n_r_old, 2, -1 + r_loc = old_radius(ir) + r_above = old_radius(ir-1) + If (abs(r_above - r_loc)/r_loc .lt. tol) Then ! the current r_loc is about to be repeated + ncheby_old(idom) = (n_r_old - ir + 1) - nr_in_lower_domains + nr_in_lower_domains = nr_in_lower_domains + ncheby_old(idom) + idom = idom + 1 + Endif + Enddo + + ! For the final (top) domain, the formula is easy + ! Note that for only one domain, this logic should still work + ! (For one domain, idom will still = 1 and nr_in_lower_domains still = 0 + ncheby_old(idom) = n_r_old - nr_in_lower_domains + End Subroutine Get_Old_Ncheby End Module Checkpointing From b013a4e33981bcdde95de387d1d18cc80ec2b21a Mon Sep 17 00:00:00 2001 From: Loren Date: Thu, 20 Jun 2024 12:07:49 -0700 Subject: [PATCH 7/8] added exit routine in the single-domain degradation case --- src/Physics/Checkpointing.F90 | 229 +++++++++++++++++++++++----------- 1 file changed, 158 insertions(+), 71 deletions(-) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index db06b6c0..cab7df5a 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -248,7 +248,7 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Integer :: read_magnetism = 0, read_hydro = 0 Integer, Allocatable :: rinds(:), gpars(:,:), read_var(:) Real*8 :: dt_pars(3),dt,new_dt - Real*8, Allocatable :: old_radius(:), radius_old_loc(:) + Real*8, Allocatable :: old_radius(:), radius_old(:), radius_old_loc(:) Real*8, Allocatable :: tempfield1(:,:,:,:), tempfield2(:,:,:,:) Character*120 :: autostring, dstring, iterstring, access_type Character*256 :: grid_file, checkfile, input_file @@ -582,109 +582,70 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif - ! Loop over the domains to set the Chebyshev coefficients, - ! possibly interpolating in radius for each subdomain - ! The fields are easy, since they are stored in spectral (Chebyshev) space - ! The AB terms are a bit more complicated, since they are currently stored in physical (radius) space - ! When we change the checkpointing format, should also store AB terms in cheby-space - irmin_old = n_r_old - irmax_old = n_r_old - ncheby_old(1) + 1 - irmin = n_r - irmax = n_r - ncheby(1) + 1 - Do idom = 1, ndomains + If (ndomains .eq. 1) Then ! Keep everything the same as it was before + ! interpolation in multiple domains was enabled + ! NOW, if n_r_old and grid_type_old are the same, we can copy chtkmp%p1b into abterms and ! fields. Otherwise, we need to interpolate onto the current grid - ! And by "interpolate," we mean set the lowest-order Chebyshev coefficients of - ! the new checkpoint using the old checkpoint, and zero out the higher-order coefficients. - ! For abterms, we will transform to spectral space, set the lowest-order coefficients and - ! and transform back - - ! Loop from inner domain to outer domain - n_r_old_loc = ncheby_old(idom) - n_r_loc = ncheby(idom) - - ! Decrement the rmin/rmax indices (which go down with increasing radius/idom - ! for each subdomain after the first - If (idom .gt. 1) Then - irmin_old = irmin_old - ncheby_old(idom - 1) - irmax_old = irmax_old - ncheby_old(idom) - irmin = irmin - ncheby(idom - 1) - irmax = irmax - ncheby(idom) - Endif - - ! Interpolate if necessary (in each domain separately) - If ((n_r_old_loc .ne. n_r_loc) .or. (grid_type_old .ne. grid_type) ) Then + ! When we change the checkpointing format, should also store AB terms in cheby-space + If ((n_r_old .ne. n_r) .or. (grid_type_old .ne. grid_type) ) Then ! Interpolate ! We will assume the user kept the same radial domain bounds. ! If they have not, this will end badly. If (my_rank .eq. 0) Then Call stdout%print(' ') - If (ndomains .gt. 1) Then ! if more than one domain, specify - ! which one we are in - Write(szstr,'(i13)')idom - Call stdout%print('------ In domain: '//TRIM(szstr)) - Endif Call stdout%print('------ Radial grid has changed.') Call stdout%print('------ Interpolating onto new grid.') Write(szstr,'(i13)')grid_type_old Call stdout%print('------ Old grid_type: '//TRIM(szstr)) Write(szstr,'(i13)')grid_type Call stdout%print('------ Current grid_type: '//TRIM(szstr)) - Write(szstr,'(i13)')n_r_old_loc + Write(szstr,'(i13)')n_r_old Call stdout%print('------ Old N_R: '//TRIM(szstr)) - Write(szstr,'(i13)')n_r_loc + Write(szstr,'(i13)')n_r Call stdout%print('------ Current N_R: '//TRIM(szstr)) Call stdout%print(' ') Endif - If (n_r_old_loc .lt. n_r_loc) Then + If (n_r_old .lt. n_r) Then + ! The fields are OK - they are already in chebyshev space - ! Just set the lowest-order coefficients from the checkpoint - fields(irmax:irmax+n_r_old_loc-1,:,:,1:numfields) = chktmp%p1b(irmax_old:irmin_old,:,:,1:numfields) + fields(:,:,:,1:numfields) = chktmp%p1b(:,:,:,1:numfields) ! The AB terms are stored in physical space (in radius). - ! They need to be transformed, lowest-order coefficients copied, and transformed back. + ! They need to be transformed, coefficients copied, and transformed back.. ! First, we need to initialize the old chebyshev grid. - Allocate(radius_old_loc(1:n_r_old_loc)) - Call cheby_info%Init(radius_old_loc,radius(irmin),radius(irmax)) ! We assume that the domain bounds do not change + Allocate(radius_old(1:n_r_old)) + Call cheby_info%Init(radius_old,rmin,rmax) ! We assume that rmax and rmin do not change fcount(:,:) = numfields - If (first_time_interpolating) Then ! Can only initialize chktmp2 once - Call chktmp2%init(field_count = fcount, config = 'p1a') - first_time_interpolating = .False. - Endif + Call chktmp2%init(field_count = fcount, config = 'p1a') Call chktmp2%construct('p1a') chktmp2%p1a(:,:,:,:) = 0.0d0 ! Allocate tempfield1, tempfield2 lb = lbound(chktmp%p1b,3) ub = ubound(chktmp%p1b,3) - Allocate(tempfield1(1:n_r_old_loc,1:2,lb:ub,1)) - Allocate(tempfield2(1:n_r_old_loc,1:2,lb:ub,1)) - - ! For each field, put the physical AB terms from the checkpoint in tempfield1 - ! use "tospec4d" to put the spectral tempfield1 coefficients in tempfield2 - ! copy the tempfield2 spectral coefficients into chktmp2 - ! (thus interpolating, since chktmp2 will have zeros for the higher-order coefficients) - Do i = 1, numfields + Allocate(tempfield1(1:n_r_old,1:2,lb:ub,1)) + Allocate(tempfield2(1:n_r_old,1:2,lb:ub,1)) + + Do i = 1, numfields tempfield1(:,:,:,:) = 0.0d0 tempfield2(:,:,:,:) = 0.0d0 - tempfield1(:,:,:,1) = chktmp%p1b(irmax_old:irmin_old,:,:,numfields+i) - Call cheby_info%tospec4d(tempfield1,tempfield2) - chktmp2%p1a(irmax:irmax+n_r_old_loc-1,:,:,i) = tempfield2(:,:,:,1) + tempfield1(1:n_r_old,:,:,1) = chktmp%p1b(1:n_r_old,:,:,numfields+i) + call cheby_info%tospec4d(tempfield1,tempfield2) + chktmp2%p1a(1:n_r_old,:,:,i) = tempfield2(1:n_r_old,:,:,1) Enddo DeAllocate(tempfield1,tempfield2) - ! chktmp2 now has the interpolated abterms from the checkpoint in spectral space - ! use "From_Spectral" to transform these fields into physical space and copy into abterms + Call chktmp2%construct('p1b') !Normal transform(p1a,p1b) Call gridcp%From_Spectral(chktmp2%p1a,chktmp2%p1b) - abterms(irmax:irmin,:,:,1:numfields) = chktmp2%p1b(irmax:irmin,:,:,1:numfields) - ! clean up + abterms(:,:,:,1:numfields) = chktmp2%p1b(:,:,:,1:numfields) Call cheby_info%destroy() Call chktmp2%deconstruct('p1a') Call chktmp2%deconstruct('p1b') - Deallocate(radius_old_loc) + Deallocate(radius_old) Else ! Rayleigh doesn't currently support degrading radial resolution--exit now If (my_rank .eq. 0) Then Call stdout%print('ERROR: Rayleigh currently does not support degrading radial resolution.') @@ -694,16 +655,142 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Call pfi%exit() Stop Endif - Else ! n_r_old_loc = n_r_loc - ! no interpolation is needed, we just copy the checkpoint into fields and abterms - fields(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,1:numfields) - abterms(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,numfields+1:numfields*2) + + Else + + ! Interpolation is complete, now we just copy into the other arrays + fields(:,:,:,1:numfields) = chktmp%p1b(:,:,:,1:numfields) + abterms(:,:,:,1:numfields) = chktmp%p1b(:,:,:,numfields+1:numfields*2) Endif - Enddo + Call chktmp%deconstruct('p1b') + DeAllocate(old_radius) + + Else ! ndomains > 1: we need to loop over domains and (maybe) interpolate + + ! Loop over the domains to set the Chebyshev coefficients, + ! possibly interpolating in radius for each subdomain + ! The fields are easy, since they are stored in spectral (Chebyshev) space + ! The AB terms are a bit more complicated, since they are currently stored in physical (radius) space + ! When we change the checkpointing format, should also store AB terms in cheby-space + irmin_old = n_r_old + irmax_old = n_r_old - ncheby_old(1) + 1 + irmin = n_r + irmax = n_r - ncheby(1) + 1 + Do idom = 1, ndomains + ! NOW, if n_r_old and grid_type_old are the same, we can copy chtkmp%p1b into abterms and + ! fields. Otherwise, we need to interpolate onto the current grid + ! And by "interpolate," we mean set the lowest-order Chebyshev coefficients of + ! the new checkpoint using the old checkpoint, and zero out the higher-order coefficients. + ! For abterms, we will transform to spectral space, set the lowest-order coefficients and + ! and transform back + + ! Loop from inner domain to outer domain + n_r_old_loc = ncheby_old(idom) + n_r_loc = ncheby(idom) + + ! Decrement the rmin/rmax indices (which go down with increasing radius/idom + ! for each subdomain after the first + If (idom .gt. 1) Then + irmin_old = irmin_old - ncheby_old(idom - 1) + irmax_old = irmax_old - ncheby_old(idom) + irmin = irmin - ncheby(idom - 1) + irmax = irmax - ncheby(idom) + Endif + + ! Interpolate if necessary (in each domain separately) + If ((n_r_old_loc .ne. n_r_loc) .or. (grid_type_old .ne. grid_type) ) Then + ! Interpolate + ! We will assume the user kept the same radial domain bounds. + ! If they have not, this will end badly. + If (my_rank .eq. 0) Then + Call stdout%print(' ') + If (ndomains .gt. 1) Then ! if more than one domain, specify + ! which one we are in + Write(szstr,'(i13)')idom + Call stdout%print('------ In domain: '//TRIM(szstr)) + Endif + Call stdout%print('------ Radial grid has changed.') + Call stdout%print('------ Interpolating onto new grid.') + Write(szstr,'(i13)')grid_type_old + Call stdout%print('------ Old grid_type: '//TRIM(szstr)) + Write(szstr,'(i13)')grid_type + Call stdout%print('------ Current grid_type: '//TRIM(szstr)) + Write(szstr,'(i13)')n_r_old_loc + Call stdout%print('------ Old N_R: '//TRIM(szstr)) + Write(szstr,'(i13)')n_r_loc + Call stdout%print('------ Current N_R: '//TRIM(szstr)) + Call stdout%print(' ') + Endif + + If (n_r_old_loc .lt. n_r_loc) Then + ! The fields are OK - they are already in chebyshev space + ! Just set the lowest-order coefficients from the checkpoint + fields(irmax:irmax+n_r_old_loc-1,:,:,1:numfields) = chktmp%p1b(irmax_old:irmin_old,:,:,1:numfields) + + ! The AB terms are stored in physical space (in radius). + ! They need to be transformed, lowest-order coefficients copied, and transformed back. + ! First, we need to initialize the old chebyshev grid. + Allocate(radius_old_loc(1:n_r_old_loc)) + Call cheby_info%Init(radius_old_loc,radius(irmin),radius(irmax)) ! We assume that the domain bounds do not change + fcount(:,:) = numfields + If (first_time_interpolating) Then ! Can only initialize chktmp2 once + Call chktmp2%init(field_count = fcount, config = 'p1a') + first_time_interpolating = .False. + Endif + Call chktmp2%construct('p1a') + chktmp2%p1a(:,:,:,:) = 0.0d0 + ! Allocate tempfield1, tempfield2 + lb = lbound(chktmp%p1b,3) + ub = ubound(chktmp%p1b,3) + Allocate(tempfield1(1:n_r_old_loc,1:2,lb:ub,1)) + Allocate(tempfield2(1:n_r_old_loc,1:2,lb:ub,1)) + + ! For each field, put the physical AB terms from the checkpoint in tempfield1 + ! use "tospec4d" to put the spectral tempfield1 coefficients in tempfield2 + ! copy the tempfield2 spectral coefficients into chktmp2 + ! (thus interpolating, since chktmp2 will have zeros for the higher-order coefficients) + Do i = 1, numfields + tempfield1(:,:,:,:) = 0.0d0 + tempfield2(:,:,:,:) = 0.0d0 + tempfield1(:,:,:,1) = chktmp%p1b(irmax_old:irmin_old,:,:,numfields+i) + Call cheby_info%tospec4d(tempfield1,tempfield2) + chktmp2%p1a(irmax:irmax+n_r_old_loc-1,:,:,i) = tempfield2(:,:,:,1) + Enddo + DeAllocate(tempfield1,tempfield2) + + ! chktmp2 now has the interpolated abterms from the checkpoint in spectral space + ! use "From_Spectral" to transform these fields into physical space and copy into abterms + Call chktmp2%construct('p1b') + !Normal transform(p1a,p1b) + Call gridcp%From_Spectral(chktmp2%p1a,chktmp2%p1b) + abterms(irmax:irmin,:,:,1:numfields) = chktmp2%p1b(irmax:irmin,:,:,1:numfields) + + ! clean up + Call cheby_info%destroy() + Call chktmp2%deconstruct('p1a') + Call chktmp2%deconstruct('p1b') + Deallocate(radius_old_loc) + Else ! Rayleigh doesn't currently support degrading radial resolution--exit now + If (my_rank .eq. 0) Then + Call stdout%print('ERROR: Rayleigh currently does not support degrading radial resolution.') + Call stdout%print('Now exiting') + Call stdout%partial_flush() + Endif + Call pfi%exit() + Stop + Endif + Else ! n_r_old_loc = n_r_loc + ! no interpolation is needed, we just copy the checkpoint into fields and abterms + fields(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,1:numfields) + abterms(irmax:irmin,:,:,1:numfields) = chktmp%p1b(irmax_old:irmax_old+n_r_loc-1,:,:,numfields+1:numfields*2) + + Endif + Enddo - DeAllocate(old_radius) - Call chktmp%deconstruct('p1b') + DeAllocate(old_radius) + Call chktmp%deconstruct('p1b') + Endif ! This ends the new block of code to deal with multiple domain interpolation End Subroutine Read_Checkpoint From c413b4a298688d81ef46a7ccd063a9ce602d14b4 Mon Sep 17 00:00:00 2001 From: Loren Date: Thu, 20 Jun 2024 16:26:54 -0700 Subject: [PATCH 8/8] updated changelog with multi domain interpolation note --- CHANGELOG.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d9f6c0a8..0eb8609a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm - The streamfunction equations solved by Rayleigh have been added to the documentation. \[Tami Rogers; 6-20-2024; [#533](https://github.com/geodynamics/Rayleigh/pull/533)\] +- Rayleigh now supports increasing radial resolution even when there are multiple Chebyshev domains. \[Loren Matilsky; 6-20-2024; [#537](https://github.com/geodynamics/Rayleigh/pull/537)\] + ### Changed - The logic in the configure script that determines the library-link flags for compilation on different systems has been simplified. \[Rene Gassmoeller; 6-19-2024; [#534](https://github.com/geodynamics/Rayleigh/pull/534)\]