diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 47e90727..12f96ce1 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -52,12 +52,32 @@ jobs: make fdeps && git diff --exit-code make -j make install - - name: Test + - name: Build Rayleigh v 1.3 + # Build 1.3 for testing checkpoint compatibility + run: | + # Build Rayleigh v 1.3 + git clone -b rayleigh-1.3.0 https://github.com/geodynamics/Rayleigh.git + mv Rayleigh Rayleigh_1.3 + cd Rayleigh_1.3 + ./configure -debian-mkl ${{ matrix.mpi }} --FFLAGS_DBG='-O0 -g -fbounds-check -ffpe-trap=invalid,zero,overflow -fbacktrace -ffixed-line-length-132 -Wall' + + #make fdeps && git diff --exit-code + make -j + make install + - name: Checkpoint Test + run: | + # Checkpoint Test + cd "$GITHUB_WORKSPACE" + cd ./tests/checkpoint + sh ./run_test.sh + + - name: Chebyshev Test run: | #Primary test (Chebyshev Mode) cd "$GITHUB_WORKSPACE"/tests/c2001_case0 mpirun -np 4 ../../bin/rayleigh.dbg - + - name: Finite-difference Test + run: | # Finite-difference test (uniform grid) cd "$GITHUB_WORKSPACE"/tests/c2001_case0_FD_uniform mpirun -np 4 ../../bin/rayleigh.dbg @@ -69,27 +89,33 @@ jobs: # J2011 steady mhd test # cd "$GITHUB_WORKSPACE"/tests/j2011_steady_mhd_minimal # mpirun -np 4 ../../bin/rayleigh.dbg - + - name: Generic Input Test + run: | # Generic input test cd "$GITHUB_WORKSPACE" sh ./tests/generic_input/run_test.sh - + - name: Chi Scalar Test + run: | # chi scalar test cd "$GITHUB_WORKSPACE" sh ./tests/chi_scalar/run_test.sh - + - name: Custom Reference Test + run: | # custom reference state test cd "$GITHUB_WORKSPACE" sh ./tests/custom_reference/run_test.sh - + - name: Chi Custom Reference Test + run: | # chi custom reference state test cd "$GITHUB_WORKSPACE" sh ./tests/chi_custom_reference/run_test.sh - + - name: Coupled BC Test + run: | # coupled bc test cd "$GITHUB_WORKSPACE" sh ./tests/coupled_bcs/run_test.sh - + - name: Vforce Diagnostics Test + run: | # vforce diagnostics regression test cd "$GITHUB_WORKSPACE" sh ./tests/vforce_diagnostics/run_test.sh diff --git a/src/IO/Parallel_IO.F90 b/src/IO/Parallel_IO.F90 index ace556fc..d1c7ebb7 100644 --- a/src/IO/Parallel_IO.F90 +++ b/src/IO/Parallel_IO.F90 @@ -180,7 +180,7 @@ Subroutine Initialize_IO_Buffer(self,grid_pars, nvals, mpi_tag, & averaging_weights, nrec, skip, & write_timestamp, averaging_axes, & spectral, mode, l_values, cache_spectral, & - spec_comp, lmax_in, full_cache) + spec_comp, nr_in, lmax_in, full_cache) Implicit None Class(io_buffer) :: self Integer, Intent(In) :: grid_pars(1:,1:) @@ -191,6 +191,7 @@ Subroutine Initialize_IO_Buffer(self,grid_pars, nvals, mpi_tag, & Logical, Intent(In), Optional :: full_cache, spec_comp Integer, Intent(In), Optional :: averaging_axes(3) Integer, Intent(In), Optional :: mode + Integer, Intent(In), Optional :: nr_in Integer, Intent(In), Optional :: lmax_in Integer :: adim(2), in_lmax Real*8, Allocatable :: avg_weights(:,:) diff --git a/src/Physics/Checkpointing.F90 b/src/Physics/Checkpointing.F90 index 46b85a68..2c8e71d4 100755 --- a/src/Physics/Checkpointing.F90 +++ b/src/Physics/Checkpointing.F90 @@ -37,10 +37,10 @@ Module Checkpointing ! Simple Checkpointing Module ! Uses MPI-IO to split writing of files amongst rank zero processes from each row Implicit None - Type(SphericalBuffer) :: chktmp, chktmp2, bctmp + Type(SphericalBuffer) :: chktmp, chktmp2, bctmp,abterms_cheby Integer, private :: numfields Integer, private :: check_err_off = 100 ! Checkpoint errors report in range 100-200. - Integer, private :: checkpoint_version = 2 + Integer, private :: checkpoint_version = 3 Integer,private :: checkpoint_tag = 425 Character*120 :: checkpoint_prefix ='nothing' Character*6 :: auto_fmt = '(i2.2)' ! Format code for quicksaves @@ -151,6 +151,11 @@ Subroutine Initialize_Checkpointing() Call checkpoint_buffer%Init(gpars, mpi_tag=checkpoint_tag, & spectral=.true., cache_spectral = .true., spec_comp = .true.) DeAllocate(gpars) + + ! Buffer to hold Adams-Bashforth terms in Chebyshev format + nfs(:) = numfields + if (magnetism) nfs(:) = numfields+1 ! magnetic mode has one extra field space used to combined derivatives + Call abterms_cheby%init(field_count = nfs, config = 'p1a') End Subroutine Initialize_Checkpointing Subroutine Write_Checkpoint(abterms,iteration,dt,new_dt,elapsed_time, input_file) @@ -165,10 +170,23 @@ Subroutine Write_Checkpoint(abterms,iteration,dt,new_dt,elapsed_time, input_file endian_tag=314 Call chktmp%construct('p1a') chktmp%config = 'p1a' - !Copy the RHS into chtkmp - Call Get_All_RHS(chktmp%p1a) - chktmp%p1a(:,:,:,numfields+1:numfields*2) = abterms(:,:,:,1:numfields) - !Now we want to move from p1a to s2a (rlm space) + + !Copy the RHS (contains state variables at current timestep) into chtkmp + Call Get_All_RHS(chktmp%p1a) + + + If ( (chebyshev) .and. (ndomains .eq. 1) ) Then + ! Convert the AB terms into chebyshev space and copy them into the buffer + ! For now, if ndomains > 1, we do not perform this conversion + Call abterms_cheby%construct('p1a') + abterms_cheby%config='p1a' + !chktmp%p1a(:,:,:,numfields+1:numfields*2) = abterms(:,:,:,1:numfields) ! earlier version -- physical space storage + Call gridcp%to_spectral(abterms,abterms_cheby%p1a) + chktmp%p1a(:,:,:,numfields+1:numfields*2) = abterms_cheby%p1a(:,:,:,1:numfields) + Call abterms_cheby%deconstruct('p1a') + Endif + + !Move checkpoint buffer from p1a to s2a (rlm space) Call chktmp%reform() If (ItIsTimeForAQuickSave) Then @@ -188,7 +206,6 @@ Subroutine Write_Checkpoint(abterms,iteration,dt,new_dt,elapsed_time, input_file checkfile = Trim(my_path)//trim(checkpoint_prefix)//'/'//trim(checkpoint_suffix(i)) Call checkpoint_buffer%cache_data_spectral(chktmp%s2a,i) Call checkpoint_buffer%write_data(filename=checkfile, clear_existing = .true.) - Enddo Call chktmp%deconstruct('s2a') @@ -262,7 +279,7 @@ 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 + nsubmax), fcount(3,2), version + Integer :: old_pars(8 + 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 @@ -374,6 +391,7 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) If (legacy_format .or. (endian_tag .eq. 314)) Then If (.not. legacy_format) Then Read(funit)version + Write(6,*)'Current Checkpoint Version: ', version Read(funit)n_r_old Else Read(funit)n_r_old @@ -441,7 +459,8 @@ 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) + old_pars(8) = version + old_pars(9:9+nsubmax-1) = ncheby_old(1:nsubmax) dt_pars(1) = dt dt_pars(2) = new_dt dt_pars(3) = checkpoint_time @@ -458,13 +477,14 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) l_max_old = old_pars(3) checkpoint_iter = old_pars(4) last_auto = old_pars(5) + version = old_pars(8) If (old_pars(6) .eq. 2) Then under_slash='_' legacy_format=.true. Endif - ncheby_old(1:nsubmax) = old_pars(8:8+nsubmax-1) + ncheby_old(1:nsubmax) = old_pars(9:9+nsubmax-1) If (old_pars(7) .ne. 0) Then ! Something is wrong with this checkpoint. @@ -566,6 +586,10 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) DeAllocate(rinds) DeAllocate(gpars) Endif + + ! Load the state vector fields and the AB terms. + ! During the read-in process below, old fields are stored in + ! indices 1:n_r_old of first dimension of chktmp%p1b. Do i = 1, numfields*2 If (read_var(i) .eq. 1) Then checkfile = trim(checkpoint_prefix)//under_slash//trim(checkpoint_suffix(i)) @@ -573,7 +597,6 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Call checkpoint_inbuffer%grab_data_spectral(chktmp%s2b,i) Endif Enddo - Call chktmp%reform() ! move to p1b If (.not. legacy_format) Then @@ -601,92 +624,131 @@ Subroutine Read_Checkpoint(fields, abterms,iteration,read_pars) Endif - 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 - ! 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(' ') - Endif + If ( (version .ge. 3) .and. chebyshev ) Then + ! For version 3+, the AB terms are stored in Chebyshev space + Call abterms_cheby%construct('p1b') + abterms_cheby%config='p1b' + abterms_cheby%p1b = 0d0 + ub = Min(n_r,n_r_old) + abterms_cheby%p1b(1:ub,:,:,1:numfields)=chktmp%p1b(1:ub,:,:,numfields+1:numfields*2) + 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) + If (ndomains .eq. 1) Then + ! Regardless of the checkpoint version or resolution, we can load the state vector + ! fields directly from the checkpoint buffer because that portion of the buffer is + ! already in the expected format. + ub = Min(n_r,n_r_old) + fields = 0d0 + fields(1:ub,:,:,1:numfields) = chktmp%p1b(1:ub,:,:,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. + If (.not. chebyshev) Then + abterms(:,:,:,1:numfields) = chktmp%p1b(:,:,:,numfields+1:numfields*2) + + Else If (version .ge. 3) Then + ! For checkpoint version 3 and later, the AB terms are also stored in + ! Chebyshev space, but they need to be in grid space. + ! Even if the resolution changes, no further work needs to be done. + !Need to convert from Chebyshev space to physical space + Call gridcp%From_Spectral(abterms_cheby%p1b,abterms) + Call abterms_cheby%deconstruct('p1b') + Else + ! For checkpoint versions 1 and 2, AB terms were stored in physical space. + ! If the grid resolution changed, we need to perform a Chebyshev interpolation. + ! If not, we can load directly from the checkpoint buffer. + ! Due to the storage convention for versions 1 and 2, interpolation + ! can only be carried out if n_r_old < n_r. + If (n_r_old .eq. n_r) Then + abterms(:,:,:,1:numfields) = chktmp%p1b(:,:,:,numfields+1:numfields*2) + Else If (n_r_old .lt. n_r) Then + Write(6,*)'Inside interpolation branch' + ! (1) 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 + Call cheby_info%Init(radius_old,rmin,rmax) ! rmax and rmin cannot change + + ! (2) Create some buffer space to store the old radial profiles and + ! the new Chebyshev coefficients. 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) + lb = lbound(chktmp%p1b,3) ! describes bounds of lm mode index 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)) + ! (3) Convert to the new Chebyshev grid Do i = 1, numfields tempfield1(:,:,:,:) = 0.0d0 tempfield2(:,:,:,:) = 0.0d0 + ! Subsample smaller-grid fields from larger-grid buffer tempfield1(1:n_r_old,:,:,1) = chktmp%p1b(1:n_r_old,:,:,numfields+i) - call cheby_info%tospec4d(tempfield1,tempfield2) + ! Transform + Call cheby_info%tospec4d(tempfield1,tempfield2) + + ! Grab Chebyshev coefficients up to old n_max chktmp2%p1a(1:n_r_old,:,:,i) = tempfield2(1:n_r_old,:,:,1) Enddo - DeAllocate(tempfield1,tempfield2) + Call gridcp%From_Spectral(chktmp2%p1a,abterms) - - Call chktmp2%construct('p1b') - !Normal transform(p1a,p1b) - Call gridcp%From_Spectral(chktmp2%p1a,chktmp2%p1b) - - abterms(:,:,:,1:numfields) = chktmp2%p1b(:,:,:,1:numfields) + ! clean up + DeAllocate(tempfield1,tempfield2) Call cheby_info%destroy() - Call chktmp2%deconstruct('p1a') - Call chktmp2%deconstruct('p1b') Deallocate(radius_old) - Else ! Rayleigh doesn't currently support degrading radial resolution--exit now + Call chktmp2%deconstruct('p1a') + Else 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%print(' ') + Call stdout%print('******************************** ERROR *********************************') + Call stdout%print('------ Rayleigh does not support degrading radial resolution') + Call stdout%print('------ for checkpoints generated using version 1.3.0 or earlier.') + Call stdout%print(' ') + Call stdout%print('------ In order to degrade the radial resolution of this model:') + Call stdout%print('------ (1) Restart from the current checkpoint using this version of') + Call stdout%print('------ Rayleigh and maintain the current radial resolution.') + Call stdout%print('------ (2) Generate a new checkpoint at the current resolution.') + Call stdout%print('------ This checkpoint will be saved using the updated format.') + Call stdout%print('------ (3) Restart using the new checkpoint and desired radial resolution.') + Call stdout%print(' ') Call stdout%partial_flush() Endif Call pfi%exit() Stop - Endif + Endif ! (n_r_old .eq. n_r) - Else + Endif ! (version .ge. 3) - ! 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) + ! clean up + Call chktmp%deconstruct('p1b') + Deallocate(old_radius) + ! Finally, print a message if we interpolated + If ( (my_rank .eq. 0) .and. (n_r .ne. n_r_old) ) 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(' ') + If (n_r_old .gt. n_r) Then + Call stdout%print('------ ******************* WARNING *****************************') + Call stdout%print('------ Boundary conditions are not formally satisfied by the') + Call stdout%print('------ checkpoint save state when N_R is decreased upon restart.') + Call stdout%print('------ Drastically degrading the radial resolution may lead') + Call stdout%print('------ to unexpected results.') + Call stdout%print(' ') + Endif Endif - Call chktmp%deconstruct('p1b') - DeAllocate(old_radius) Else ! ndomains > 1: we need to loop over domains and (maybe) interpolate - + ! Note that ndomains > 1 is not possible when running in finite-difference mode. ! 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 diff --git a/src/Physics/Fields.F90 b/src/Physics/Fields.F90 index 72c99c97..730a7f59 100755 --- a/src/Physics/Fields.F90 +++ b/src/Physics/Fields.F90 @@ -452,17 +452,17 @@ Subroutine Initialize_Field_Structure() ! These following code should pretty much never be modified by the user. if (.not. magnetism) then - wsfcount(1,2) = n_equations + n_active_scalars + n_passive_scalars - wsfcount(2,2) = n_equations + n_active_scalars + n_passive_scalars - wsfcount(3,2) = n_equations + n_active_scalars + n_passive_scalars + wsfcount(1,2) = n_equations + wsfcount(2,2) = n_equations + wsfcount(3,2) = n_equations else emfr = avar emftheta = cvar emfphi = avar+1 ! seven RHS's (plus scalars) go back for the solve (1 field is differentiated and combined at the end) - wsfcount(1,2) = n_equations + 1 + n_active_scalars + n_passive_scalars - wsfcount(2,2) = n_equations + 1 + n_active_scalars + n_passive_scalars - wsfcount(3,2) = n_equations + 1 + n_active_scalars + n_passive_scalars + wsfcount(1,2) = n_equations + 1 + wsfcount(2,2) = n_equations + 1 + wsfcount(3,2) = n_equations + 1 endif diff --git a/tests/checkpoint/compare_check.py b/tests/checkpoint/compare_check.py new file mode 100644 index 00000000..a9d995ae --- /dev/null +++ b/tests/checkpoint/compare_check.py @@ -0,0 +1,58 @@ +from rayleigh_diagnostics import G_Avgs, SPH_Modes +import numpy as np + +def chisq(a1,a2): + b = (a1-a2)**2 + c = np.sum(b)/np.sum(a1*a1) + return c**0.5 + +dirs = ['old_format', 'new_format', 'new_format_same_res', 'new_format_up_res', 'new_format_down_res'] +iters = ['00000100', '00000150', '00000200','00000250'] +tols = [1e-10, 1e-10, 1e-2, 1e-2] +ndirs = len(dirs) +for i in range(1,ndirs): + #print('=========================') + chisqs = [] + old = dirs[i-1] + new = dirs[i] + istring = iters[i-1] + #print(old,new,istring) + + sphm_o = SPH_Modes(istring,path=old+'/SPH_Modes/') + sphm_n = SPH_Modes(istring,path=new+'/SPH_Modes/') + ga_o = G_Avgs(istring,path=old+'/G_Avgs/') + ga_n = G_Avgs(istring,path=new+'/G_Avgs/') + + mx_ref = np.max([sphm_o.vals.real**2,sphm_o.vals.imag**2])**0.5 + for j, l in enumerate(sphm_o.lvals): + for m in range(0,l+1): + val_o = sphm_o.vals[m,j,0,0,:] + val_n = sphm_n.vals[m,j,0,0,:] + c = chisq(val_o.real,val_n.real) + + mx = np.max([val_o.real**2,val_n.imag**2])**0.5 + chisqs.append(c*mx/mx_ref) + + if (m > 0): # m = 0 has no imaginary component + c = chisq(val_o.imag,val_n.imag) + chisqs.append(c*mx/mx_ref) + if ((l == 4000) and (m ==4)): + fig,ax = plt.subplots(ncols=2,figsize=(10,5)) + ax[0].plot(val_o.real) + ax[0].plot(val_n.real) + ax[1].plot(val_o.imag) + ax[1].plot(val_n.imag) + plt.show() + sph_csq = np.max(chisqs) + + val_o = ga_o.vals[:,0] + val_n = ga_n.vals[:,0] + ga_csq = chisq(val_o,val_n) + + #print('sph modes: ', sph_csq) + #print('G_Avgs: ', ga_csq) + if ( (ga_csq > tols[i-1]) or (sph_csq > tols[i-1]) ): + print('Checkpoint Test Error: Time series do not agree.') + exit(1) +print('Checkpoint Test Passed') +exit(0) diff --git a/tests/checkpoint/old_format/main_input b/tests/checkpoint/old_format/main_input new file mode 100644 index 00000000..3065b2ee --- /dev/null +++ b/tests/checkpoint/old_format/main_input @@ -0,0 +1,77 @@ +&problemsize_namelist + n_r=48 + n_theta = 64 + nprow = 2 + npcol = 2 + aspect_ratio = 0.35d0 + shell_depth = 1.0d0 +/ +&numerical_controls_namelist +/ +&physical_controls_namelist +! benchmark_mode = 1 +! benchmark_integration_interval = 100 +! benchmark_report_interval = 5000 + rotation = .True. + magnetism = .false. + viscous_heating = .false. + ohmic_heating = .false. + advect_reference_state = .false. +/ +&temporal_controls_namelist + max_time_step = 1.0d-4 + max_iterations = 200 + checkpoint_interval = 25 + cflmin = 0.4d0 + cflmax = 0.6d0 +/ +&io_controls_namelist +/ +&output_namelist + +! New equatorial slices output + +sph_mode_ell = 0,4 +sph_mode_levels = 1 ! outer boundary +sph_mode_values = 507 ! radial entropy gradient +sph_mode_frequency = 10 +sph_mode_nrec = 5 + +globalavg_values = 401 ! kinetic energy +globalavg_frequency = 10 +globalavg_nrec = 5 + + +full3d_values = 64 ! temperature +full3d_frequency = 9000000 +/ + +&Boundary_Conditions_Namelist +no_slip_boundaries = .true. +strict_L_Conservation = .false. +dtdr_bottom = 0.0d0 +T_Top = 0.0d0 +T_Bottom = 1.0d0 +fix_tvar_top = .true. +fix_tvar_bottom = .true. +/ +&Initial_Conditions_Namelist +init_type=1 ! Benchmark init +temp_amp = 1.0d-3 +temp_w = 0.01d4 +restart_iter=-1 +conductive_profile=.true. +/ +&Test_Namelist +/ +&Reference_Namelist +Ekman_Number = 1.0d-3 +Rayleigh_Number = 1.0d5 +Prandtl_Number = 1.0d0 +Magnetic_Prandtl_Number = 5.0d0 +reference_type = 1 +heating_type = 0 ! No heating +gravity_power = 1.0d0 ! g ~ radius +/ +&Transport_Namelist +/ diff --git a/tests/checkpoint/run_test.sh b/tests/checkpoint/run_test.sh new file mode 100644 index 00000000..8626bbde --- /dev/null +++ b/tests/checkpoint/run_test.sh @@ -0,0 +1,93 @@ +#!/usr/bin/env bash + +export RA_ROOT=../../.. +export rayleigh_v13=$RA_ROOT/Rayleigh_1.3/bin/rayleigh.dbg +export rayleigh_dbg=$RA_ROOT/bin/rayleigh.dbg + +#cd tests/checkpoint +cp ../../post_processing/rayleigh_diagnostics.py . +# Generate a checkpoint in the old format using Rayleigh v1.3 +cd old_format +mpirun -np 4 $rayleigh_v13 -niter 100 +cd .. + +# Restart from the old checkpoint format at increased resolution using Rayleigh 1.3. +lastdir=old_format +newdir=old_format_up + +mkdir $newdir +mkdir $newdir/Checkpoints +cp $lastdir/main_input $newdir/. +cp -r $lastdir/Checkpoints/00000050 $newdir/Checkpoints/. + +cd $newdir +sed -i 's/n_r=48/n_r=64/g' main_input +sed -i 's/init_type=7/init_type=-1/g' main_input +sed -i 's/restart_iter=-1/restart_iter=50/g' main_input +mpirun -np 4 $rayleigh_v13 -niter 100 +cd .. + + +# Restart from the old checkpoint format at same resolution using current Rayleigh. +lastdir=old_format +newdir=new_format + +mkdir $newdir +mkdir $newdir/Checkpoints +cp $lastdir/main_input $newdir/. +cp -r $lastdir/Checkpoints/00000050 $newdir/Checkpoints/. + +cd $newdir +#sed -i 's/n_r=48/n_r=64/g' main_input +sed -i 's/init_type=7/init_type=-1/g' main_input +sed -i 's/restart_iter=-1/restart_iter=50/g' main_input +mpirun -np 4 $rayleigh_dbg -niter 100 +cd .. + + +# Restart from the new checkpoint format at same resolution. +lastdir=new_format +newdir=new_format_same_res +mkdir $newdir +mkdir $newdir/Checkpoints +cp $lastdir/main_input $newdir/. +cp -r $lastdir/Checkpoints/00000100 $newdir/Checkpoints/. + +cd $newdir +#sed -i 's/n_r=48/n_r=64/g' main_input +#sed -i 's/init_type=1/init_type=-1/g' main_input +sed -i 's/restart_iter=50/restart_iter=100/g' main_input +mpirun -np 4 $rayleigh_dbg -niter 100 +cd .. + +# Restart from the new checkpoint format at increased resolution. +lastdir=new_format_same_res +newdir=new_format_up_res +mkdir $newdir +mkdir $newdir/Checkpoints +cp $lastdir/main_input $newdir/. +cp -r $lastdir/Checkpoints/00000150 $newdir/Checkpoints/. + +cd $newdir +sed -i 's/n_r=48/n_r=64/g' main_input +#sed -i 's/init_type=1/init_type=-1/g' main_input +sed -i 's/restart_iter=100/restart_iter=150/g' main_input +mpirun -np 4 $rayleigh_dbg -niter 100 +cd .. + +# Restart from the new checkpoint format and degrade resolution. +lastdir=new_format_up_res +newdir=new_format_down_res +mkdir $newdir +mkdir $newdir/Checkpoints +cp $lastdir/main_input $newdir/. +cp -r $lastdir/Checkpoints/00000200 $newdir/Checkpoints/. + +cd $newdir +sed -i 's/n_r=64/n_r=48/g' main_input +#sed -i 's/init_type=1/init_type=-1/g' main_input +sed -i 's/restart_iter=150/restart_iter=200/g' main_input +mpirun -np 4 $rayleigh_dbg -niter 100 +cd .. + +python3 compare_check.py