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
25 procedure:: pre_update
33 SUBROUTINE init (self, link)
35 class(
link_t),
pointer :: link
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
45 integer,
save :: axis = 3
46 namelist /gravity_params/ on, constant, axis, position, mass, minrad
48 call trace%begin (
'selfgravity_t%init')
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
59 self%patch => task2patch(link%task)
61 self%const_grav = constant
63 self%position = position
66 if (constant /= 0.0)
then 67 self%grav => const_grav
68 else if (mass > -1.0)
then 71 call io%abort(
"gravity_mod%init:: neither constant nor mass is set.")
73 if (.not.
allocated(self%patch%force_per_unit_mass))
then 75 allocate (self%patch%force_per_unit_mass(gn(1),gn(2),gn(3),3))
76 self%patch%force_per_unit_mass = 0.0
84 SUBROUTINE pre_update(self)
88 END SUBROUTINE pre_update
92 SUBROUTINE const_grav(self)
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
101 SUBROUTINE var_grav(self)
104 integer:: ix, iy, iz, i
106 real(8):: x, y, z, xs, ys, zs, xp, yp, zp, ds2(3)
107 real(kind=KindScalarVar):: newton
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)
113 ds2(i) = m(i)%h(self%patch%idx%px+i-1) *m(i)%d
115 do iz=m(3)%lb,m(3)%ub
117 zp = z - self%position(3)
119 do iy=m(2)%lb,m(2)%ub
121 yp = y - self%position(2)
123 do ix=m(1)%lb,m(1)%ub
125 xp = x - self%position(1)
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
148 END SUBROUTINE var_grav
Define code units, in terms of (here) CGS units.
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Add gravity to any solver.
Fundamental constants in CGS and SI units.
Module with list handling for generic class task_t objects.