Merge pull request #488 from bhishekmanek/NewInitialConditions

Initial conditions using checkpoint files added to user supplied files

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

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
@@ -126,11 +125,26 @@ 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
+
+        ! 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)
+            If (init_type .eq. -2) Then
+                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
+                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
         Endif
 
-
         !////////////////////////////////////
         ! Initialize the hydro variables
 
@@ -224,8 +238,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 +636,34 @@ 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
+        call tempfield%init(field_count = fcount, config = 'p1b')
+        call tempfield%construct('p1b')
+        call tempfield2%init(field_count = fcount, config = 'p1b')
+        call tempfield2%construct('p1b')
+
+        call get_rhs(field_index,tempfield2%p1b(:,:,:,1))
+
+        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)