From d8a44a006f74f5ea3185506a3370876161f6ac38 Mon Sep 17 00:00:00 2001 From: Bhishek Manek Date: Tue, 3 Oct 2023 11:43:25 -0600 Subject: [PATCH 1/3] Updating Initial Conditions This file includes init_type = -2 for adding entropy perturbations from checkpoint and an user input file. Also included is a magnetic_init_type = -2 that adds checkpoints to user input files of poloidal and toroidal component of magnetic fields. --- src/Physics/Initial_Conditions.F90 | 65 +++++++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/src/Physics/Initial_Conditions.F90 b/src/Physics/Initial_Conditions.F90 index cbd896f6..8c9b339b 100755 --- a/src/Physics/Initial_Conditions.F90 +++ b/src/Physics/Initial_Conditions.F90 @@ -126,10 +126,38 @@ Subroutine Initialize_Fields() Call stdout%print(" ---- Magnetic Init Type : RESTART ") Endif Endif - If ( (init_type .eq. -1) .or. ( magnetism .and. (magnetic_init_type .eq. -1) ) ) Then + + If ( (init_type .lt. 0) .or. ( magnetism .and. (magnetic_init_type .lt. 0) ) ) Then Call restart_from_checkpoint(restart_iter) + print*,'Init type < 0; Restarting from Checkpoint' + ! Bhishek + If (init_type .eq. -2) Then + ! Add entropy perturbation + print*,'Adding Entropy perturbation file' + Call add_to_field(tvar, t_init_file) + Endif + If (magnetic_init_type .eq. -2) Then + print*,'Adding poloidal and toroidal componets to a checkpoint' + Call add_to_field(cvar, c_init_file) + Call add_to_field(avar, a_init_file) + Endif + ! Bhishek + Endif + + !NICK + If (init_type .eq. -2) Then + print*, 'Init type -2; Restart with checkpoint' + !Call restart_from_checkpoint(restart_iter) + ! Add entropy perturbation + print*,'Init type -2; Adding Entropy perturbation file' + Call add_to_field(tvar, t_init_file) Endif + If (magnetic_init_type .eq. -2) Then + ! Add A and C perturbations + Call add_to_field(cvar, c_init_file) + Call add_to_field(avar, a_init_file) + Endif !//////////////////////////////////// ! Initialize the hydro variables @@ -224,8 +252,10 @@ Subroutine Restart_From_Checkpoint(iteration) If (magnetism) Then + If (init_type .eq. -2) rpars(1) = 1 If (init_type .eq. -1) rpars(1) = 1 If (magnetic_init_type .eq. -1) rpars(2) = 1 + If (magnetic_init_type .eq. -2) rpars(2) = 1 !If both variable types are not read in, an euler_step is taken on restart prod = rpars(1)*rpars(2) @@ -620,6 +650,39 @@ subroutine magnetic_file_init() end subroutine magnetic_file_init + !NICK + subroutine add_to_field(field_index, field_file) + ! initialize magnetic variables from generic input files + Implicit None + Integer, Intent(In) :: field_index + Character*120, Intent(In) :: field_file + Integer :: fcount(3,2) + Type(SphericalBuffer) :: tempfield, tempfield2 + fcount(:,:) = 1 + print*,'Success 1' + call tempfield%init(field_count = fcount, config = 'p1b') + print*,'Success 2' + call tempfield%construct('p1b') + print*,'Success 3' + call tempfield2%init(field_count = fcount, config = 'p1b') + call tempfield2%construct('p1b') + print*,'Success 4' + + call get_rhs(field_index,tempfield2%p1b(:,:,:,1)) + print*,'Success 5' + + if (trim(field_file) .ne. '__nothing__') then + call read_input(field_file, 1, tempfield) + + tempfield%p1b = tempfield%p1b+tempfield2%p1b + call set_rhs(field_index, tempfield%p1b(:,:,:,1)) + end if + + + call tempfield%deconstruct('p1b') + call tempfield2%deconstruct('p1b') + + end subroutine add_to_field !////////////////////////////////////////////////////////////////////////////////// ! Diffusion Init (for linear solve development) From 65e5496d815bb153c4ca56e1d9f5b2f17854b4d5 Mon Sep 17 00:00:00 2001 From: Bhishek Manek Date: Thu, 5 Oct 2023 18:19:02 -0700 Subject: [PATCH 2/3] Adding Set_RHS subroutine to Linear_Solve --- src/Math_Layer/Linear_Solve.F90 | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/Math_Layer/Linear_Solve.F90 b/src/Math_Layer/Linear_Solve.F90 index 70ad739b..33b96e65 100755 --- a/src/Math_Layer/Linear_Solve.F90 +++ b/src/Math_Layer/Linear_Solve.F90 @@ -641,6 +641,29 @@ Subroutine Set_RHS(eqid,set_to) Endif End Subroutine Set_RHS + Subroutine Get_RHS(eqid,get_from) + ! Set the RHS of equation eqid to the value of set_to + Implicit None + Real*8, Intent(InOut) :: get_from(:,:,:) + Integer, Intent(In) :: eqid + Integer :: istart,iend, ind + + + If (n_modes .gt. 0) Then + !Primary equation object always has the full, allocated rhs. + ind = eqid + if (.not. equation_set(1,eqid)%primary) then + ind = equation_set(1,eqid)%links(1) + endif + ! Individual RHS's inhabit row ranges defined by rowblock and ndim1 + istart = equation_set(1,eqid)%rowblock+1 + iend = istart+ndim1-1 + + + get_from(1:ndim1,:,:)=equation_set(1,ind)%rhs(istart:iend,:,:) + Endif + End Subroutine Get_RHS + Subroutine Zero_RHS(eqid) ! Set the RHS of equation eqid to zero ! This is useful for when we want to remove explicity time-dependence From 29e58a0af51cd1e5120357dfb2f86c8f706103bf Mon Sep 17 00:00:00 2001 From: Bhishek Manek Date: Thu, 5 Oct 2023 22:29:11 -0700 Subject: [PATCH 3/3] Cleaing up Initial Conditions file --- src/Physics/Initial_Conditions.F90 | 35 +++++++----------------------- 1 file changed, 8 insertions(+), 27 deletions(-) diff --git a/src/Physics/Initial_Conditions.F90 b/src/Physics/Initial_Conditions.F90 index 8c9b339b..e641723c 100755 --- a/src/Physics/Initial_Conditions.F90 +++ b/src/Physics/Initial_Conditions.F90 @@ -113,7 +113,6 @@ Subroutine Initialize_Fields() ! The equation set RHS's stays allocated throughout - it is effectively how we save the AB terms. Call Allocate_RHS(zero_rhs=.true.) - !//////////////////////////////////////// ! Read in checkpoint files as appropriate If (init_type .eq. -1) Then @@ -127,36 +126,23 @@ Subroutine Initialize_Fields() Endif Endif + ! Take care of IC's that require reading from a checkpoint (if init_type and magnetic_init_type = -1) + ! If (init_type or/and magnetic_init_type = -2), read the checkpoint and add a user defined field to it for IC If ( (init_type .lt. 0) .or. ( magnetism .and. (magnetic_init_type .lt. 0) ) ) Then Call restart_from_checkpoint(restart_iter) - print*,'Init type < 0; Restarting from Checkpoint' - ! Bhishek If (init_type .eq. -2) Then - ! Add entropy perturbation - print*,'Adding Entropy perturbation file' + If (my_rank .eq. 0) Then + Call stdout%print(" ---- Hydro Init Type : ADD USER FILE TO CHECKPOINT ") + Endif Call add_to_field(tvar, t_init_file) Endif If (magnetic_init_type .eq. -2) Then - print*,'Adding poloidal and toroidal componets to a checkpoint' + If (my_rank .eq. 0) Then + Call stdout%print(" ---- Magnetic Init Type : ADD USER FILE TO CHECKPOINT ") + Endif Call add_to_field(cvar, c_init_file) Call add_to_field(avar, a_init_file) Endif - ! Bhishek - Endif - - !NICK - If (init_type .eq. -2) Then - print*, 'Init type -2; Restart with checkpoint' - !Call restart_from_checkpoint(restart_iter) - ! Add entropy perturbation - print*,'Init type -2; Adding Entropy perturbation file' - Call add_to_field(tvar, t_init_file) - Endif - - If (magnetic_init_type .eq. -2) Then - ! Add A and C perturbations - Call add_to_field(cvar, c_init_file) - Call add_to_field(avar, a_init_file) Endif !//////////////////////////////////// @@ -659,17 +645,12 @@ subroutine add_to_field(field_index, field_file) Integer :: fcount(3,2) Type(SphericalBuffer) :: tempfield, tempfield2 fcount(:,:) = 1 - print*,'Success 1' call tempfield%init(field_count = fcount, config = 'p1b') - print*,'Success 2' call tempfield%construct('p1b') - print*,'Success 3' call tempfield2%init(field_count = fcount, config = 'p1b') call tempfield2%construct('p1b') - print*,'Success 4' call get_rhs(field_index,tempfield2%p1b(:,:,:,1)) - print*,'Success 5' if (trim(field_file) .ne. '__nothing__') then call read_input(field_file, 1, tempfield)