<commit_before>module grav_sources_module

  implicit none

  private

  public add_grav_source

contains

! :::
! ::: ------------------------------------------------------------------
! :::

    subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, &
                               uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, &
                               grav, gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, &
                               lo,hi,dt,dx,E_added)

      use bl_constants_module, only: M_PI
      use meth_params_module, only : NVAR, URHO, UMX, UMY, UMZ, UEDEN, grav_source_type, rot_freq
      use prob_params_module, only: coord_type
      use probdata_module, only: center

      implicit none

      integer lo(3), hi(3)
      integer uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3
      integer  uout_l1, uout_l2, uout_l3, uout_h1, uout_h2, uout_h3
      integer  gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3

      double precision  uin( uin_l1: uin_h1, uin_l2: uin_h2, uin_l3: uin_h3,NVAR)
      double precision uout(uout_l1:uout_h1,uout_l2:uout_h2,uout_l3:uout_h3,NVAR)
      double precision grav(  gv_l1:  gv_h1,  gv_l2:  gv_h2,  gv_l3:  gv_h3,3)
      double precision dt, dx(3)
      double precision E_added

      double precision :: rho
      double precision :: SrU, SrV, SrW, SrE
      double precision :: rhoInv
      double precision :: old_rhoeint, new_rhoeint, old_ke, new_ke, old_re
      double precision :: x, y, z
      integer          :: i, j, k

      double precision, parameter :: TWO_PI = 2.d0 * M_PI
      double precision :: omega

      ! should come in through a module
      double precision :: problo(3)


      if (coord_type == 0) then
         omega = TWO_PI*rot_freq
      else
         call bl_error("Error:: Rotate_3d.f90 :: unknown coord_type")
      endif

      problo(:) = 0.0d0

      ! hack: do both gravity and rotation here, so we can couple them.
      ! we need to make sure that Rotate_3d.f90 does nothing.

      do k = lo(3),hi(3)
         z = problo(3) + dx(3)*(float(k)+0.5d0) - center(3)

         do j = lo(2),hi(2)
            y = problo(2) + dx(2)*(float(j)+0.5d0) - center(2)

            do i = lo(1),hi(1)
               x = problo(1) + dx(1)*(float(i)+0.5d0) - center(1)


               rho    = uin(i,j,k,URHO)
               rhoInv = 1.0d0 / rho

               ! momenta

               ! gravity
               SrU = rho * grav(i,j,k,1)
               SrV = rho * grav(i,j,k,2)
               SrW = rho * grav(i,j,k,3)

               ! rotation
               SrU = SrU + 2.d0*uin(i,j,k,UMY)*omega + rho*omega**2*x
               SrV = SrV - 2.d0*uin(i,j,k,UMX)*omega + rho*omega**2*y

               ! update the momenta
               uout(i,j,k,UMX)   = uout(i,j,k,UMX) + SrU * dt
               uout(i,j,k,UMY)   = uout(i,j,k,UMY) + SrV * dt
               uout(i,j,k,UMZ)   = uout(i,j,k,UMZ) + SrW * dt

               
               ! energy

               ! Src = rho u dot g, evaluated with all quantities at t^n
               SrE = uin(i,j,k,UMX) * grav(i,j,k,1) + &
                     uin(i,j,k,UMY) * grav(i,j,k,2) + &
                     uin(i,j,k,UMZ) * grav(i,j,k,3)

               ! rotation
               SrE = SrE + omega**2*(uin(i,j,k,UMX)*x + uin(i,j,k,UMY)*y)

               uout(i,j,k,UEDEN) = uout(i,j,k,UEDEN) + SrE * dt

            enddo
         enddo
      enddo

      end subroutine add_grav_source

end module grav_sources_module

! :::
! ::: ------------------------------------------------------------------
! :::

      subroutine ca_corrgsrc(lo,hi, &
                             gold,gold_l1,gold_l2,gold_l3,gold_h1,gold_h2,gold_h3, &
                             gnew,gnew_l1,gnew_l2,gnew_l3,gnew_h1,gnew_h2,gnew_h3, &
                             uold,uold_l1,uold_l2,uold_l3,uold_h1,uold_h2,uold_h3, &
                             unew,unew_l1,unew_l2,unew_l3,unew_h1,unew_h2,unew_h3, &
                             dt,dx,E_added)

      use bl_constants_module, only: M_PI
      use meth_params_module, only : NVAR, URHO, UMX, UMY, UMZ, UEDEN, grav_source_type, rot_freq
      use prob_params_module, only: coord_type
      use probdata_module, only: center

      implicit none

      integer lo(3),hi(3)
      integer gold_l1,gold_l2,gold_l3,gold_h1,gold_h2,gold_h3
      integer gnew_l1,gnew_l2,gnew_l3,gnew_h1,gnew_h2,gnew_h3
      integer uold_l1,uold_l2,uold_l3,uold_h1,uold_h2,uold_h3
      integer unew_l1,unew_l2,unew_l3,unew_h1,unew_h2,unew_h3
      double precision   gold(gold_l1:gold_h1,gold_l2:gold_h2,gold_l3:gold_h3,3)
      double precision   gnew(gnew_l1:gnew_h1,gnew_l2:gnew_h2,gnew_l3:gnew_h3,3)
      double precision  uold(uold_l1:uold_h1,uold_l2:uold_h2,uold_l3:uold_h3,NVAR)
      double precision  unew(unew_l1:unew_h1,unew_l2:unew_h2,unew_l3:unew_h3,NVAR)
      double precision  dt,dx(3),E_added
      double precision :: x, y, z
      integer i,j,k

      double precision Sx, Sy, Sz
      double precision SrEcorr
      double precision rhoo
      double precision rhon

      double precision rhooinv, rhoninv
      double precision old_ke, old_rhoeint, old_re
      double precision new_ke, new_rhoeint

      double precision, parameter :: TWO_PI = 2.d0 * M_PI
      double precision :: omega

      ! should come in through a module
      double precision :: problo(3)


      if (coord_type == 0) then
         omega = TWO_PI*rot_freq
      else
         call bl_error("Error:: Rotate_3d.f90 :: unknown coord_type")
      endif

      problo(:) = 0.0d0

      do k = lo(3),hi(3)
         z = problo(3) + dx(3)*(float(k)+0.5d0) - center(3)

         do j = lo(2),hi(2)
            y = problo(2) + dx(2)*(float(j)+0.5d0) - center(2)

            do i = lo(1),hi(1)
               x = problo(1) + dx(1)*(float(i)+0.5d0) - center(1)


               rhoo    = uold(i,j,k,URHO)
               rhooinv = 1.0d0 / uold(i,j,k,URHO)

               rhon    = unew(i,j,k,URHO)
               rhoninv = 1.0d0 / unew(i,j,k,URHO)


               ! momenta source corrections
               Sx = unew(i,j,k,UMX) + &
                    0.5d0*dt*(rhon*gnew(i,j,k,1) - rhoo*gold(i,j,k,1)) + &
                    0.5d0*dt*(rhon*omega**2*x - rhoo*omega**2*x) - dt*uold(i,j,k,UMY)*omega

               Sy = unew(i,j,k,UMY) + &
                    0.5d0*dt*(rhon*gnew(i,j,k,2) - rhoo*gold(i,j,k,2)) + &
                    0.5d0*dt*(rhon*omega**2*y - rhoo*omega**2*y) + dt*uold(i,j,k,UMX)*omega

               Sz = unew(i,j,k,UMZ) + &<commit_msg>only deal with rotation here in the correction.<commit_after>module grav_sources_module

  implicit none

  private

  public add_grav_source

contains

! :::
! ::: ------------------------------------------------------------------
! :::

    subroutine add_grav_source(uin,uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3, &
                               uout,uout_l1,uout_l2,uout_l3,uout_h1,uout_h2,uout_h3, &
                               grav, gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3, &
                               lo,hi,dt,dx,E_added)

      use bl_constants_module, only: M_PI
      use meth_params_module, only : NVAR, URHO, UMX, UMY, UMZ, UEDEN, grav_source_type, rot_freq
      use prob_params_module, only: coord_type
      use probdata_module, only: center

      implicit none

      integer lo(3), hi(3)
      integer uin_l1,uin_l2,uin_l3,uin_h1,uin_h2,uin_h3
      integer  uout_l1, uout_l2, uout_l3, uout_h1, uout_h2, uout_h3
      integer  gv_l1, gv_l2, gv_l3, gv_h1, gv_h2, gv_h3

      double precision  uin( uin_l1: uin_h1, uin_l2: uin_h2, uin_l3: uin_h3,NVAR)
      double precision uout(uout_l1:uout_h1,uout_l2:uout_h2,uout_l3:uout_h3,NVAR)
      double precision grav(  gv_l1:  gv_h1,  gv_l2:  gv_h2,  gv_l3:  gv_h3,3)
      double precision dt, dx(3)
      double precision E_added

      double precision :: rho
      double precision :: SrU, SrV, SrW, SrE
      double precision :: rhoInv
      double precision :: old_rhoeint, new_rhoeint, old_ke, new_ke, old_re
      double precision :: x, y, z
      integer          :: i, j, k

      double precision, parameter :: TWO_PI = 2.d0 * M_PI
      double precision :: omega

      ! should come in through a module
      double precision :: problo(3)


      if (coord_type == 0) then
         omega = TWO_PI*rot_freq
      else
         call bl_error("Error:: Rotate_3d.f90 :: unknown coord_type")
      endif

      problo(:) = 0.0d0

      ! hack: do both gravity and rotation here, so we can couple them.
      ! we need to make sure that Rotate_3d.f90 does nothing.

      do k = lo(3),hi(3)
         z = problo(3) + dx(3)*(float(k)+0.5d0) - center(3)

         do j = lo(2),hi(2)
            y = problo(2) + dx(2)*(float(j)+0.5d0) - center(2)

            do i = lo(1),hi(1)
               x = problo(1) + dx(1)*(float(i)+0.5d0) - center(1)


               rho    = uin(i,j,k,URHO)
               rhoInv = 1.0d0 / rho

               ! momenta

               ! gravity
               SrU = rho * grav(i,j,k,1)
               SrV = rho * grav(i,j,k,2)
               SrW = rho * grav(i,j,k,3)

               ! rotation
               !SrU = SrU + 2.d0*uin(i,j,k,UMY)*omega + rho*omega**2*x
               !SrV = SrV - 2.d0*uin(i,j,k,UMX)*omega + rho*omega**2*y

               ! update the momenta
               uout(i,j,k,UMX)   = uout(i,j,k,UMX) + SrU * dt
               uout(i,j,k,UMY)   = uout(i,j,k,UMY) + SrV * dt
               uout(i,j,k,UMZ)   = uout(i,j,k,UMZ) + SrW * dt

               
               ! energy

               ! Src = rho u dot g, evaluated with all quantities at t^n
               SrE = uin(i,j,k,UMX) * grav(i,j,k,1) + &
                     uin(i,j,k,UMY) * grav(i,j,k,2) + &
                     uin(i,j,k,UMZ) * grav(i,j,k,3)

               ! rotation
               !SrE = SrE + omega**2*(uin(i,j,k,UMX)*x + uin(i,j,k,UMY)*y)

               uout(i,j,k,UEDEN) = uout(i,j,k,UEDEN) + SrE * dt

            enddo
         enddo
      enddo

      end subroutine add_grav_source

end module grav_sources_module

! :::
! ::: ------------------------------------------------------------------
! :::

      subroutine ca_corrgsrc(lo,hi, &
                             gold,gold_l1,gold_l2,gold_l3,gold_h1,gold_h2,gold_h3, &
                             gnew,gnew_l1,gnew_l2,gnew_l3,gnew_h1,gnew_h2,gnew_h3, &
                             uold,uold_l1,uold_l2,uold_l3,uold_h1,uold_h2,uold_h3, &
                             unew,unew_l1,unew_l2,unew_l3,unew_h1,unew_h2,unew_h3, &
                             dt,dx,E_added)

      use bl_constants_module, only: M_PI
      use meth_params_module, only : NVAR, URHO, UMX, UMY, UMZ, UEDEN, grav_source_type, rot_freq
      use prob_params_module, only: coord_type
      use probdata_module, only: center

      implicit none

      integer lo(3),hi(3)
      integer gold_l1,gold_l2,gold_l3,gold_h1,gold_h2,gold_h3
      integer gnew_l1,gnew_l2,gnew_l3,gnew_h1,gnew_h2,gnew_h3
      integer uold_l1,uold_l2,uold_l3,uold_h1,uold_h2,uold_h3
      integer unew_l1,unew_l2,unew_l3,unew_h1,unew_h2,unew_h3
      double precision   gold(gold_l1:gold_h1,gold_l2:gold_h2,gold_l3:gold_h3,3)
      double precision   gnew(gnew_l1:gnew_h1,gnew_l2:gnew_h2,gnew_l3:gnew_h3,3)
      double precision  uold(uold_l1:uold_h1,uold_l2:uold_h2,uold_l3:uold_h3,NVAR)
      double precision  unew(unew_l1:unew_h1,unew_l2:unew_h2,unew_l3:unew_h3,NVAR)
      double precision  dt,dx(3),E_added
      double precision :: x, y, z
      integer i,j,k

      double precision Sx, Sy, Sz
      double precision SrEcorr
      double precision rhoo
      double precision rhon

      double precision rhooinv, rhoninv
      double precision old_ke, old_rhoeint, old_re
      double precision new_ke, new_rhoeint

      double precision, parameter :: TWO_PI = 2.d0 * M_PI
      double precision :: omega

      ! should come in through a module
      double precision :: problo(3)


      if (coord_type == 0) then
         omega = TWO_PI*rot_freq
      else
         call bl_error("Error:: Rotate_3d.f90 :: unknown coord_type")
      endif

      problo(:) = 0.0d0

      print *, 'in corrgsrc, rot_freq = ', rot_freq

      do k = lo(3),hi(3)
         z = problo(3) + dx(3)*(float(k)+0.5d0) - center(3)

         do j = lo(2),hi(2)
            y = problo(2) + dx(2)*(float(j)+0.5d0) - center(2)

            do i = lo(1),hi(1)
               x = problo(1) + dx(1)*(float(i)+0.5d0) - center(1)


               rhoo    = uold(i,j,k,URHO)
               rhooinv = 1.0d0 / uold(i,j,k,URHO)

               rhon    = unew(i,j,k,URHO)
               rhoninv = 1.0d0 / unew(i,j,k,URHO)


               ! momenta source corrections
               Sx = unew(i,j,k,UMX) + &
                    0.5d0*dt*(rhon*gnew(i,j,k,1) - rhoo*gold(i,j,k,1)) + &
                    0.5d0*dt*(rhon*omega**2*x - rhoo*omega**2*x) - dt*uold(i,j,k,UMY)*omega

               Sy = unew(i,j,k,UMY) + &
                    0.5d0*dt*(rhon*gnew(i,j,k,2) - rhoo*gold(i,j,k,2)) + &
                    0.5d0*dt*(rhon*omega**2*y - rhoo*omega**2*y) + dt*uold(i,j,k,UMX)*omega

               Sz = unew(i,j,k,UMZ) + &
{'fraction_non_alphanumeric': 0.14143183549124144, 'fraction_numerical': 0.03777608530083778, 'mean_word_length': 1.918, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6326743', 'n_tokens_mistral': 5604, 'n_tokens_neox': 5369, 'n_words': 902}