DISPATCH
gravity_mod.f90
1 !===============================================================================
2 !> Add gravity to any solver
3 !===============================================================================
4 MODULE gravity_mod
5  USE patch_mod
6  USE units_mod
7  USE scaling_mod
8  USE link_mod
9  USE io_mod
10  USE kinds_mod
11  USE trace_mod
12  implicit none
13  private
14  type, public:: gravity_t
15  procedure(const_grav), pointer :: grav
16  real(kind=KindScalarVar) :: const_grav = 0.0
17  real(kind=KindScalarVar) :: mass = -1.0
18  real(8) :: minrad = 0.0
19  real(8) :: position(3) = 0.0 ! centre of gravity
20  integer :: axis = 3
21  class(patch_t), pointer:: patch
22  logical:: on=.false.
23  contains
24  procedure:: init
25  procedure:: pre_update
26  end type
27  type(gravity_t), public:: gravity
28 CONTAINS
29 
30 !===============================================================================
31 !> Initialize gravity
32 !===============================================================================
33 SUBROUTINE init (self, link)
34  class(gravity_t):: self
35  class(link_t), pointer :: link
36  ! ............................................................................
37  class(patch_t), pointer:: patch
38  integer :: iostat, gn(3)
39  real(kind=KindScalarVar), save :: constant = 0.0
40  logical, save :: first_time=.true.
41  logical, save :: on = .false.
42  real(kind=KindScalarVar), save:: mass = -1.0
43  real(8), save:: minrad = 0.0
44  real(8), save :: position(3) = 0.0 ! centre of gravity
45  integer, save :: axis = 3
46  namelist /gravity_params/ on, constant, axis, position, mass, minrad
47  !-----------------------------------------------------------------------------
48  call trace%begin ('selfgravity_t%init')
49  !$omp critical (input_cr)
50  if (first_time) then
51  first_time = .false.
52  rewind(io%input)
53  read (io%input, gravity_params, iostat=iostat)
54  if (io%master) write (io%output, gravity_params)
55  if (constant/=0.0) constant = constant * scaling%t**2 / scaling%l
56  end if
57  !$omp end critical (input_cr)
58  if (on) then
59  self%patch => task2patch(link%task)
60  self%on = on
61  self%const_grav = constant
62  self%axis = axis
63  self%position = position
64  self%mass = mass
65  self%minrad = minrad
66  if (constant /= 0.0) then
67  self%grav => const_grav
68  else if (mass > -1.0) then
69  self%grav => var_grav
70  else
71  call io%abort("gravity_mod%init:: neither constant nor mass is set.")
72  end if
73  if (.not.allocated(self%patch%force_per_unit_mass)) then
74  gn = self%patch%gn
75  allocate (self%patch%force_per_unit_mass(gn(1),gn(2),gn(3),3))
76  self%patch%force_per_unit_mass = 0.0
77  end if
78  end if
79  call trace%end()
80 END SUBROUTINE init
81 
82 !===============================================================================
83 !===============================================================================
84 SUBROUTINE pre_update(self)
85  class(gravity_t) :: self
86  ! ----------------------------------------------------------------------------
87  call self%grav()
88 END SUBROUTINE pre_update
89 
90 !===============================================================================
91 !===============================================================================
92 SUBROUTINE const_grav(self)
93  class(gravity_t) :: self
94  ! ----------------------------------------------------------------------------
95  self%patch%force_per_unit_mass(:,:,:,self%axis) = &
96  self%patch%force_per_unit_mass(:,:,:,self%axis) + self%const_grav
97 END SUBROUTINE const_grav
98 
99 !===============================================================================
100 !===============================================================================
101 SUBROUTINE var_grav(self)
102  class(gravity_t) :: self
103  ! ----------------------------------------------------------------------------
104  integer:: ix, iy, iz, i
105  real(8):: rp
106  real(8):: x, y, z, xs, ys, zs, xp, yp, zp, ds2(3)
107  real(kind=KindScalarVar):: newton
108  ! ............................................................................
109  newton = self%mass*scaling%grav
110  associate(m=>self%patch%mesh)
111  associate(r1=>m(1)%r, r2=>m(2)%r, r3=>m(3)%r)
112  do i=1,3
113  ds2(i) = m(i)%h(self%patch%idx%px+i-1) *m(i)%d
114  end do
115  do iz=m(3)%lb,m(3)%ub
116  z = m(3)%p + r3(iz)
117  zp = z - self%position(3)
118  zs = zp + ds2(3)
119  do iy=m(2)%lb,m(2)%ub
120  y = m(2)%p + r2(iy)
121  yp = y - self%position(2)
122  ys = yp + ds2(2)
123  do ix=m(1)%lb,m(1)%ub
124  x = m(1)%p + r1(ix)
125  xp = x - self%position(1)
126  xs = xp + ds2(1)
127  !-----------------------------------------------------------------------
128  ! The force of gravity is included explicitly, and is oriented
129  ! away from the center.
130  !-----------------------------------------------------------------------
131  rp = sqrt(xs**2+yp**2+zp**2)
132  rp = max(rp,self%minrad)
133  self%patch%force_per_unit_mass(ix,iy,iz,1) = &
134  self%patch%force_per_unit_mass(ix,iy,iz,1) - newton*xs/rp**3
135  rp = sqrt(xp**2+ys**2+zp**2)
136  rp = max(rp,self%minrad)
137  self%patch%force_per_unit_mass(ix,iy,iz,2) = &
138  self%patch%force_per_unit_mass(ix,iy,iz,2) - newton*ys/rp**3
139  rp = sqrt(xp**2+yp**2+zs**2)
140  rp = max(rp,self%minrad)
141  self%patch%force_per_unit_mass(ix,iy,iz,3) = &
142  self%patch%force_per_unit_mass(ix,iy,iz,3) - newton*zs/rp**3
143  end do
144  end do
145  end do
146  end associate
147  end associate
148 END SUBROUTINE var_grav
149 
150 END MODULE gravity_mod
Define code units, in terms of (here) CGS units.
Definition: scaling_mod.f90:4
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
Add gravity to any solver.
Definition: gravity_mod.f90:4
Fundamental constants in CGS and SI units.
Definition: units_mod.f90:4
Definition: io_mod.f90:4