gravity_pointmass.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: gravity_pointmass.f90 #
5!# #
6!# Copyright (C) 2007-2024 #
7!# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8!# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
10!# #
11!# This program is free software; you can redistribute it and/or modify #
12!# it under the terms of the GNU General Public License as published by #
13!# the Free Software Foundation; either version 2 of the License, or (at #
14!# your option) any later version. #
15!# #
16!# This program is distributed in the hope that it will be useful, but #
17!# WITHOUT ANY WARRANTY; without even the implied warranty of #
18!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
19!# NON INFRINGEMENT. See the GNU General Public License for more #
20!# details. #
21!# #
22!# You should have received a copy of the GNU General Public License #
23!# along with this program; if not, write to the Free Software #
24!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
25!# #
26!#############################################################################
38!----------------------------------------------------------------------------!
47!----------------------------------------------------------------------------!
48MODULE gravity_pointmass_mod
49 USE gravity_base_mod
50 USE boundary_base_mod
51 USE fluxes_base_mod
52 USE physics_base_mod
53 USE mesh_base_mod
54 USE logging_base_mod
55 USE marray_base_mod
56 USE marray_compound_mod
57 USE common_dict
58#ifdef PARALLEL
59#ifdef HAVE_MPI_MOD
60 USE mpi
61#endif
62#endif
63 IMPLICIT NONE
64#ifdef PARALLEL
65#ifdef HAVE_MPIF_H
66 include 'mpif.h'
67#endif
68#endif
69 !--------------------------------------------------------------------------!
70 PRIVATE
71 CHARACTER(LEN=16), PARAMETER :: potential_name(2) = (/ &
72 "Newton ", &
73 "Paczynski-Wiita " /)
74 INTEGER, PARAMETER :: newton = 1
75 INTEGER, PARAMETER :: wiita = 2
76
77 TYPE, EXTENDS(gravity_base) :: gravity_pointmass
78 CLASS(logging_base), ALLOCATABLE :: potential
79 INTEGER :: outbound
80 REAL, POINTER :: mass
81 REAL, POINTER :: accrate
82 REAL, POINTER :: massloss
83 REAL :: mdot
84 REAL :: acclimit
85 REAL :: switchon
86 REAL, DIMENSION(:,:), POINTER :: pos
87 REAL, DIMENSION(:), POINTER :: r0
88 REAL, DIMENSION(:,:,:), POINTER :: r_prim
89 REAL, DIMENSION(:,:,:,:), POINTER :: fr_prim
90 REAL, DIMENSION(:,:,:,:), POINTER :: posvec_prim
91 REAL, DIMENSION(:,:,:,:), POINTER :: posvec_prim_tmp
92 REAL, DIMENSION(:,:,:,:,:), POINTER :: fposvec_prim
93 TYPE(marray_base), ALLOCATABLE :: omega, & !< angular velocity
94 omega2
95 CONTAINS
96 PROCEDURE :: initgravity_pointmass
97 PROCEDURE :: setoutput
98 PROCEDURE :: calcpotential
99 PROCEDURE :: updategravity_single
100 PROCEDURE :: infogravity
101 PROCEDURE :: calcdiskheight_single
102 final :: finalize
103 END TYPE
104
105 !--------------------------------------------------------------------------!
106 PUBLIC :: &
107 ! types
108 gravity_pointmass, &
109 ! constants
110 newton, wiita, &
111 ! elemental functions
112 getdiskheight
113 !--------------------------------------------------------------------------!
114
115CONTAINS
116
132 SUBROUTINE initgravity_pointmass(this,Mesh,Physics,config,IO)
133 USE geometry_logcylindrical_mod, ONLY: geometry_logcylindrical
134 USE geometry_spherical_mod, ONLY: geometry_spherical
135 USE geometry_cylindrical_mod, ONLY: geometry_cylindrical
136 IMPLICIT NONE
137 !------------------------------------------------------------------------!
138 CLASS(gravity_pointmass), INTENT(INOUT) :: this
139 CLASS(mesh_base), INTENT(IN) :: Mesh
140 CLASS(physics_base), INTENT(IN) :: Physics
141 TYPE(dict_typ), POINTER :: config,IO
142 !------------------------------------------------------------------------!
143 INTEGER :: potential_number, valwrite, gtype
144 INTEGER :: err,i,j,k,l
145 REAL :: r_schwarzschild=0.0,eps
146 !------------------------------------------------------------------------!
147 CALL this%InitGravity(mesh,physics,"central point mass",config,io)
148 !\todo last dimensions are absolutely not clear if right.. What are the standing for?
149 ALLOCATE(this%potential, this%pot, this%omega, this%omega2, &
150 this%r_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
151 this%fr_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
152 this%posvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
153 this%posvec_prim_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
154 this%fposvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),&
155 ! last two entries (EAST,NORTH,TOP)x(dim1,dim2,dim3)
156 this%mass, this%accrate, this%massloss, this%pos(1,3), &
157 stat = err)
158 IF (err.NE.0) CALL this%Error("InitGravity_pointmass", "Memory allocation failed!")
159
160 ! type of the potential
161 CALL getattr(config, "potential", potential_number, newton)
162 CALL this%potential%InitLogging(potential_number,potential_name(potential_number))
163
164 ! init mesh array for potential and angular velocity
165 this%pot = marray_base(4)
166 this%omega = marray_base()
167 this%omega2 = marray_base()
168 this%pot%data1d(:) = 0.0
169 this%omega%data1d(:) = 0.0
170 this%omega2%data1d(:) = 0.0
171
172 ! set (initial) mass
173 CALL getattr(config, "mass", this%mass, 1.0)
174 CALL setattr(io, "mass", ref(this%mass))
175
176 ! cartesian position of the point mass
177 CALL getattr(config, "x", this%pos(1,1), 0.0)
178 CALL getattr(config, "y", this%pos(1,2), 0.0)
179 CALL getattr(config, "z", this%pos(1,3), 0.0)
180
181 ! Softening (e.g. for planets inside the computational domain)
182 CALL getattr(config, "softening", eps, 0.0)
183
184 ! soft switch on
185 CALL getattr(config, "switchon", this%switchon, -1.0)
186
187 ! accretion limit (e.g. eddington limit) divided by central mass
188 ! the units are therefore (e.g. SI): kg/s / central_mass = 1 / s
189 ! -1.0 == off
190 CALL getattr(config, "acclimit", this%acclimit, -1.0)
191
192 CALL setattr(io, "accrate", ref(this%accrate))
193 IF (this%acclimit.GT.0.0) &
194 CALL setattr(io, "massloss", ref(this%massloss))
195
196 SELECT CASE(potential_number)
197 CASE(newton)
198 r_schwarzschild = 0.0
199 CASE(wiita)
200 ! Schwarzschild radius
201 r_schwarzschild = 2.*physics%constants%GN * this%mass / physics%constants%C**2
202 CASE DEFAULT
203 CALL this%Error("InitGravity_pointmass", "potential must be either NEWTON or WIITA")
204 END SELECT
205
206 ! reset mass flux and massloss and register for output
207 this%mdot = 0.0
208 this%massloss = 0.0
209 this%accrate = 0.0
210
211 ! define position vector from the central object to all cell bary centers
212 ! and the corresponding distances
213 IF (all(abs(this%pos(1,:)).LE.tiny(this%pos))) THEN
214 ! no shift: point mass is located at the origin of the mesh
215 ! set position vector and inverse radius using appropriate mesh arrays
216 this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:)
217 this%r_prim(:,:,:) = mesh%radius%bcenter(:,:,:)
218 this%fr_prim(:,:,:,1) = mesh%radius%faces(:,:,:,east)
219 this%fr_prim(:,:,:,2) = mesh%radius%faces(:,:,:,north)
220 this%fr_prim(:,:,:,3) = mesh%radius%faces(:,:,:,top)
221 ELSE
222 ! shifted point mass position:
223 ! compute curvilinear components of shift vector
224 ! PLEASE NOTE: We need the shifted curvilinear components at the bary_centers and
225 ! faces in EASTERN and NORTHERN direction. This is the case because we want to
226 ! calculate the potential at these face positions further below.
227 this%posvec_prim_tmp(:,:,:,1) = this%pos(1,1)
228 this%posvec_prim_tmp(:,:,:,2) = this%pos(1,2)
229 this%posvec_prim_tmp(:,:,:,3) = this%pos(1,3)
230
231 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,east,:), &
232 this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,1,:))
233 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,north,:), &
234 this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,2,:))
235 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,top,:), &
236 this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,3,:))
237
238 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_prim_tmp,this%posvec_prim)
239
240 ! subtract the result from the position vector:
241 ! this gives you the curvilinear components of all vectors pointing
242 ! from the point mass to the bary center of any cell on the mesh
243 this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_prim(:,:,:,:)
244 this%fposvec_prim(:,:,:,1,:) = mesh%posvec%faces(:,:,:,east,:) - this%fposvec_prim(:,:,:,1,:) ! EAST
245 this%fposvec_prim(:,:,:,2,:) = mesh%posvec%faces(:,:,:,north,:) - this%fposvec_prim(:,:,:,2,:) ! NORTH
246 this%fposvec_prim(:,:,:,3,:) = mesh%posvec%faces(:,:,:,top,:) - this%fposvec_prim(:,:,:,3,:) ! TOP
247
248 this%r_prim(:,:,:) = sqrt(this%posvec_prim(:,:,:,1)**2 &
249 + this%posvec_prim(:,:,:,2)**2 + this%posvec_prim(:,:,:,3)**2)
250 this%fr_prim(:,:,:,1) = sqrt(this%fposvec_prim(:,:,:,1,1)**2 &
251 + this%fposvec_prim(:,:,:,1,2)**2 + this%fposvec_prim(:,:,:,1,3)**2) ! shifted EAST-faces
252 this%fr_prim(:,:,:,2) = sqrt(this%fposvec_prim(:,:,:,2,1)**2 &
253 + this%fposvec_prim(:,:,:,2,2)**2 + this%fposvec_prim(:,:,:,2,3)**2) ! shifted NORTH-faces
254 this%fr_prim(:,:,:,3) = sqrt(this%fposvec_prim(:,:,:,3,1)**2 &
255 +this%fposvec_prim(:,:,:,3,2)**2 + this%fposvec_prim(:,:,:,3,3)**2) ! shifted TOP-faces
256 END IF
257
258 ! Initialize gravitational acceleration and Keplerian angular velocity.
259 ! Loop over ghost cells as well to ensure that all values are well defined
260 ! and remain finite; this is of particular importance for the calculation
261 ! of the disk scale height (see below).
262!NEC$ COLLAPSE
263 DO k=mesh%KGMIN,mesh%KGMAX
264 DO j=mesh%JGMIN,mesh%JGMAX
265!NEC$ IVDEP
266 DO i=mesh%IGMIN,mesh%IGMAX
267 ! compute local Keplerian angular velocity
268 ! Since omega usually becomes infinite in the limit r -> r_prim (or a, if a > 0 )
269 ! (unless softening is enabled, i. e. eps > 0) we avoid devision by zero by first
270 ! limiting the inverse of omega to some very small value.
271 this%omega2%data3d(i,j,k) = 1./ max(10*tiny(eps), &
272 (( this%r_prim(i,j,k) * (this%r_prim(i,j,k)-r_schwarzschild)**2) + eps**3) / (physics%Constants%GN*this%mass) )
273 this%omega%data3d(i,j,k) = sqrt(this%omega2%data3d(i,j,k))
274 ! curvilinear components of the gravitational acceleration;
275 ! account for dimensionality of the acceleration vector (could be < 3)
276!NEC$ SHORTLOOP
277 DO l=1,physics%VDIM
278 this%accel%data4d(i,j,k,l) = -this%omega2%data3d(i,j,k) * this%posvec_prim(i,j,k,physics%VIDX(l))
279 END DO
280 END DO
281 END DO
282 END DO
283
284 !! \todo implement computation for pseudo-Newton Paczynski-Wiita potential
285 CALL this%CalcPotential(mesh,physics,this%mass,this%r_prim,this%fr_prim,this%pot%data4d)
286
287 CALL this%SetOutput(mesh,physics,config,io)
288
289 ! enable mass accretion by setting "outbound" to one of the four boundaries
290 ! of the computational domain (depends on mesh geometry)
291 SELECT TYPE(geo=>mesh%geometry)
292 CLASS IS(geometry_cylindrical)
293 CALL getattr(config, "outbound", this%outbound, west)
294 CLASS IS(geometry_spherical)
295 CALL getattr(config, "outbound", this%outbound, west)
296 CLASS DEFAULT
297 CALL getattr(config, "outbound", this%outbound, 0)
298 CALL this%WARNING("GravitySources_pointmass","geometry does not support accretion")
299 END SELECT
300 CALL this%InfoGravity()
301 END SUBROUTINE initgravity_pointmass
302
303 SUBROUTINE calcpotential(this,Mesh,Physics,mass,dist,dist_faces,pot)
304 IMPLICIT NONE
305 !------------------------------------------------------------------------!
306 CLASS(gravity_pointmass) :: this
307 CLASS(mesh_base), INTENT(IN) :: Mesh
308 CLASS(physics_base), INTENT(IN) :: Physics
309 REAL :: mass
310 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX), &
311 INTENT(IN) :: dist
312 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,3), &
313 INTENT(IN) :: dist_faces
314 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX,4), &
315 INTENT(OUT) :: pot
316 !------------------------------------------------------------------------!
317
318 !! \todo implement computation for pseudo-Newton Paczynski-Wiita potential
319 pot(:,:,:,1) = - physics%constants%GN * mass / dist(:,:,:)
320 pot(:,:,:,2) = - physics%constants%GN * mass / dist_faces(:,:,:,1) ! direction EAST !Mesh%radius%faces(:,:,:,EAST)
321 pot(:,:,:,3) = - physics%constants%GN * mass / dist_faces(:,:,:,2) ! direction NORTH !Mesh%radius%faces(:,:,:,NORTH)
322 pot(:,:,:,4) = - physics%constants%GN * mass / dist_faces(:,:,:,3) ! direction TOP !Mesh%radius%faces(:,:,:,TOP)
323 END SUBROUTINE
324
325 SUBROUTINE infogravity(this)
326 IMPLICIT NONE
327 !------------------------------------------------------------------------!
328 CLASS(gravity_pointmass), INTENT(IN) :: this
329 !------------------------------------------------------------------------!
330 CHARACTER(LEN=32) :: param_str
331 !------------------------------------------------------------------------!
332 WRITE (param_str,'(ES10.2)') this%mass
333 CALL this%Info(" potential: " // &
334 trim(this%potential%GetName()))
335 CALL this%Info(" initial mass: " // trim(param_str))
336 WRITE (param_str, '(3(ES10.2))') this%pos(1,1:3)
337 CALL this%Info(" cart. position: " // trim(param_str))
338 IF (this%outbound.GT.0) THEN
339 param_str = "enabled"
340 ELSE
341 param_str = "disabled"
342 END IF
343 CALL this%Info(" accretion: " // trim(param_str))
344 IF (this%outbound.GT.0.AND.this%acclimit.GT.0.0) THEN
345 WRITE (param_str,'(ES10.2)') this%acclimit
346 CALL this%Info(" accretion limit: " // trim(param_str))
347 END IF
348 IF (this%switchon.GT.0.0) THEN
349 WRITE (param_str,'(ES10.2)') this%switchon
350 CALL this%Info(" switchon time: " // trim(param_str))
351 END IF
352 END SUBROUTINE infogravity
353
354 SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
355 IMPLICIT NONE
356 !------------------------------------------------------------------------!
357 CLASS(gravity_pointmass), INTENT(INOUT) :: this
358 CLASS(mesh_base), INTENT(IN) :: Mesh
359 CLASS(physics_base), INTENT(IN) :: Physics
360 TYPE(dict_typ), POINTER :: config,IO
361 !------------------------------------------------------------------------!
362 INTEGER :: valwrite
363 !------------------------------------------------------------------------!
364 ! output cartesian positions of the point mass
365 valwrite = 0
366 CALL getattr(config, "output/position", valwrite, 0)
367 IF (valwrite .EQ. 1) &
368 CALL setattr(io, "position", this%pos)
369 valwrite = 0
370 CALL getattr(config, "output/potential", valwrite, 0)
371 IF ((valwrite .EQ. 1).AND. ALLOCATED(this%pot)) THEN
372 IF (ASSOCIATED(this%pot%data4d)) &
373 CALL setattr(io, "potential", &
374 this%pot%data4d(mesh%IMIN:mesh%IMAX,mesh%JMIN:mesh%JMAX,mesh%KMIN:mesh%KMAX,4))
375 END IF
376 END SUBROUTINE setoutput
377
383 SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
384 IMPLICIT NONE
385 !------------------------------------------------------------------------!
386 CLASS(gravity_pointmass), INTENT(INOUT) :: this
387 CLASS(mesh_base), INTENT(IN) :: Mesh
388 CLASS(physics_base), INTENT(IN) :: Physics
389 CLASS(fluxes_base), INTENT(IN) :: Fluxes
390 CLASS(marray_compound), INTENT(INOUT) :: pvar
391 REAL, INTENT(IN) :: time,dt
392 !------------------------------------------------------------------------!
393 INTEGER :: i,j,k,l
394 REAL, DIMENSION(Physics%VNUM) :: bflux
395 REAL :: scaling,oldmass,massfac,sqrmassfac,dmass,dmasslim
396 !------------------------------------------------------------------------!
397 ! update accel and omega only in case of accretion
398 IF (this%outbound.NE.0) THEN
399 ! get boundary flux
400#ifdef PARALLEL
401 bflux(:) = fluxes%GetBoundaryFlux(mesh,physics,this%outbound,mpi_comm_world)
402#else
403 bflux(:) = fluxes%GetBoundaryFlux(mesh,physics,this%outbound)
404#endif
405 ! store old mass
406 oldmass = this%mass
407 ! compute the mass flux difference
408 ! REMARK #1: bflux is the accumulated flux across the boundary
409 ! REMARK #2: mdot is NOT the accretion rate, but the
410 ! mass flux from the previous evalution (see below)
411 dmass = (bflux(physics%DENSITY) - this%mdot)
412 IF(this%acclimit.GE.0.) THEN
413 dmasslim = min(dmass,this%acclimit*dt*this%mass)
414 this%massloss = this%massloss + dmass - dmasslim
415 dmass = dmasslim
416 END IF
417 ! compute new mass and current true, i. e. possibly limited, accretion rate
418 this%mass = this%mass + dmass
419 IF (dt.GT.0.0) this%accrate = dmass/dt
420 ! scale gravitational acceleration and Keplerian angular velocity
421 ! with newmass/oldmass to account for accretion
422 massfac = this%mass/oldmass
423 sqrmassfac = sqrt(massfac)
424 this%accel%data1d(:) = this%accel%data1d(:) * massfac
425 this%pot%data1d(:) = this%pot%data1d(:) * massfac
426 this%omega%data1d(:) = this%omega%data1d(:) * sqrmassfac
427 this%mdot = bflux(physics%DENSITY)
428 END IF
429
430 ! modify acceleration during switchon phase
431 IF (time.GE.0.0.AND.time.LE.this%switchon) THEN
432 ! compute new scaling
433 scaling = sin(0.5*pi*time/this%switchon)**2
434 ! compute acceleration and scale it;
435 ! during the switchon phase accretion is disabled
436!NEC$ COLLAPSE
437 DO k=mesh%KGMIN,mesh%KGMAX
438 DO j=mesh%JGMIN,mesh%JGMAX
439!NEC$ IVDEP
440 DO i=mesh%IGMIN,mesh%IGMAX
441!NEC$ SHORTLOOP
442 DO l=1,physics%VDIM
443 this%accel%data4d(i,j,k,l) = -scaling * this%omega2%data3d(i,j,k) * this%posvec_prim(i,j,k,physics%VIDX(l))
444 END DO
445 END DO
446 END DO
447 END DO
448 END IF
449
450 END SUBROUTINE updategravity_single
451
452
456 PURE SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
457 IMPLICIT NONE
458 !------------------------------------------------------------------------!
459 CLASS(gravity_pointmass), INTENT(INOUT) :: this
460 CLASS(mesh_base), INTENT(IN) :: mesh
461 CLASS(physics_base), INTENT(IN) :: physics
462 CLASS(marray_compound), INTENT(INOUT) :: pvar
463 TYPE(marray_base), INTENT(INOUT) :: bccsound,h_ext,height
464 !------------------------------------------------------------------------!
465 ! compute disk height
466 h_ext%data1d(:) = getdiskheight(bccsound%data1d(:),this%omega%data1d(:))
467 END SUBROUTINE calcdiskheight_single
468
469 SUBROUTINE finalize(this)
470 IMPLICIT NONE
471 !------------------------------------------------------------------------!
472 TYPE(gravity_pointmass), INTENT(INOUT) :: this
473 !------------------------------------------------------------------------!
474 CHARACTER(LEN=128) :: buffer
475 !------------------------------------------------------------------------!
476 IF (this%GetRank().EQ.0) THEN
477 CALL this%Info("-------------------------------------------------------------------")
478 WRITE(buffer, "(A,(ES11.3))") " central mass: ", this%mass
479 CALL this%Info(buffer)
480 IF (this%acclimit.GT.0.0) THEN
481 WRITE(buffer, "(A,(ES11.3))") " mass loss: ", this%massloss
482 CALL this%Info(buffer)
483 END IF
484 END IF
485
486 DEALLOCATE(this%omega,this%omega2,this%r_prim,this%fr_prim, &
487 this%posvec_prim,this%posvec_prim_tmp,this%fposvec_prim,this%mass, &
488 this%accrate,this%massloss,this%pos)
489
490 IF (ALLOCATED(this%potential)) DEALLOCATE(this%potential)
491 IF (ALLOCATED(this%pot)) DEALLOCATE(this%pot)
492 CALL this%Finalize_base()
493 END SUBROUTINE finalize
494
504 ELEMENTAL FUNCTION getdiskheight(cs,omega) RESULT(height)
505 IMPLICIT NONE
506 !------------------------------------------------------------------------!
507 REAL, INTENT(IN) :: cs,omega
508 REAL :: height
509 !------------------------------------------------------------------------!
510 height = cs / omega
511 END FUNCTION getdiskheight
512
513END MODULE gravity_pointmass_mod
boundary_base_mod
Definition: boundary_base.f90:40
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
geometry_cylindrical_mod
defines properties of a 3D cylindrical mesh
Definition: geometry_cylindrical.f90:36
geometry_logcylindrical_mod
defines properties of a 3D logcylindrical mesh
Definition: geometry_logcylindrical.f90:33
geometry_spherical_mod
defines properties of a 3D spherical mesh
Definition: geometry_spherical.f90:42
gravity_base_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: gravity_base.f90:46
gravity_pointmass_mod
source terms module for gravitational acceleration due to a point mass at the center of the coordinat...
Definition: gravity_pointmass.f90:48
gravity_pointmass_mod::wiita
integer, parameter, public wiita
Definition: gravity_pointmass.f90:75
gravity_pointmass_mod::potential_name
character(len=16), dimension(2), parameter potential_name
Definition: gravity_pointmass.f90:71
gravity_pointmass_mod::updategravity_single
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
updates the gravitational acceleration of the pointmass module
Definition: gravity_pointmass.f90:384
gravity_pointmass_mod::infogravity
subroutine infogravity(this)
Definition: gravity_pointmass.f90:326
gravity_pointmass_mod::calcpotential
subroutine calcpotential(this, Mesh, Physics, mass, dist, dist_faces, pot)
Definition: gravity_pointmass.f90:304
gravity_pointmass_mod::initgravity_pointmass
subroutine initgravity_pointmass(this, Mesh, Physics, config, IO)
constructor for point mass module
Definition: gravity_pointmass.f90:133
gravity_pointmass_mod::newton
integer, parameter, public newton
Definition: gravity_pointmass.f90:74
gravity_pointmass_mod::calcdiskheight_single
pure subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
computes pressure scale height of a geometrically thin Keplerian disk
Definition: gravity_pointmass.f90:457
gravity_pointmass_mod::getdiskheight
elemental real function, public getdiskheight(cs, omega)
return pressure scale height of a geometrically thin Keplerian disk
Definition: gravity_pointmass.f90:505
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
logging_base_mod::finalize
subroutine finalize(this)
Destructor of logging_base class.
Definition: logging_base.f90:160
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
mesh_base_mod
basic mesh module
Definition: mesh_base.f90:72
mesh_base_mod::east
integer, parameter east
named constant for eastern boundary
Definition: mesh_base.f90:101
mesh_base_mod::setoutput
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
Definition: mesh_base.f90:828
mesh_base_mod::top
integer, parameter top
named constant for top boundary
Definition: mesh_base.f90:101
mesh_base_mod::north
integer, parameter north
named constant for northern boundary
Definition: mesh_base.f90:101
mesh_base_mod::west
integer, parameter west
named constant for western boundary
Definition: mesh_base.f90:101
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:432
common_dict::dict_typ
Definition: common_dict.f90:123
fluxes_base_mod::fluxes_base
Definition: fluxes_base.f90:61
geometry_cylindrical_mod::geometry_cylindrical
Definition: geometry_cylindrical.f90:41
geometry_logcylindrical_mod::geometry_logcylindrical
Definition: geometry_logcylindrical.f90:39
geometry_spherical_mod::geometry_spherical
Definition: geometry_spherical.f90:49
gravity_base_mod::gravity_base
Definition: gravity_base.f90:58
gravity_pointmass_mod::gravity_pointmass
Definition: gravity_pointmass.f90:77
logging_base_mod::logging_base
common data structure
Definition: logging_base.f90:56
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:69
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
mesh_base_mod::mesh_base
mesh data structure
Definition: mesh_base.f90:122
physics_base_mod::physics_base
Definition: physics_base.f90:69
Generated by   doxygen 1.9.4