78 REAL,
DIMENSION(:,:,:),
POINTER :: r_sec
79 REAL,
DIMENSION(:,:,:,:),
POINTER :: posvec_sec
80 REAL,
DIMENSION(:,:,:,:),
POINTER :: posvec_sec_tmp
81 REAL,
DIMENSION(:,:,:,:),
POINTER :: fr_sec
82 REAL,
DIMENSION(:,:,:,:,:),
POINTER :: fposvec_sec
83 REAL,
POINTER :: mass2
131 CALL this%InitGravity(mesh,physics,
"binary point masses",config,io)
139 CALL this%Error(
"InitGravity_binary", &
140 "Physics not supported for binary gravity term")
144 this%mass,this%massloss,this%accrate,this%mass2,this%pos(3,2),this%r0(3), &
145 this%omega2, this%omega, &
146 this%r_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
147 this%r_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
148 this%posvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
149 this%posvec_prim_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
150 this%posvec_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
151 this%posvec_sec_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
152 this%pot,this%pot_prim,this%pot_sec, &
153 this%fr_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
154 this%fr_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
155 this%fposvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),&
156 this%fposvec_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),&
159 IF (err.NE.0)
CALL this%Error(
"gravity_pointmass::InitGravity_pointmass",
"Unable allocate memory!")
167 this%pot%data1d(:) = 0.0
168 this%pot%data1d(:) = 0.0
169 this%pot_prim%data1d(:) = 0.0
170 this%pot_sec%data1d(:) = 0.0
171 this%omega2%data1d(:) = 0.0
174 WHERE (mesh%without_ghost_zones%mask1d(:))
175 this%omega%data1d(:) = 0.0
177 this%omega%data1d(:) = 1.0
182 CALL getattr(config,
"mass1", this%mass, 1.0)
183 CALL setattr(io,
"mass1", ref(this%mass))
186 CALL getattr(config,
"mass2", this%mass2, 1.0)
187 CALL setattr(io,
"mass2", ref(this%mass2))
190 CALL getattr(config,
"excentricity", this%excent, 0.0)
193 CALL getattr(config,
"semimajoraxis", this%semaaxis, -1.0)
196 CALL getattr(config,
"softening1", this%eps1, 0.0)
198 CALL getattr(config,
"softening2", this%eps2, 0.0)
201 CALL getattr(config,
"switchon1", this%switchon, -1.0)
204 CALL getattr(config,
"switchon2", this%switchon2, -1.0)
207 CALL getattr(config,
"omega_rot", this%omega_rot, 0.0)
210 CALL getattr(config,
"mesh_rot", mesh_rot, 0)
211 IF (mesh_rot.EQ.1)
THEN
212 this%omega_rot = mesh%OMEGA
213 this%r0(:) = mesh%rotcent(:)
214 IF (mesh%omega.EQ.0.0) &
215 CALL this%Warning(
"InitGravity_binary",
"Rotating mesh enabled but OMEGA set to zero in mesh")
226 this%period = 2.*pi*sqrt(this%semaaxis/(this%mass+this%mass2)/physics%constants%GN)*this%semaaxis
229 CALL this%SetOutput(mesh,physics,config,io)
230 CALL this%InfoGravity()
241 CHARACTER(LEN=32) :: mass1_str,mass2_str,excent_str,sema_str,omega_str
243 WRITE (mass1_str,
'(ES10.2)') this%mass
244 WRITE (mass2_str,
'(ES10.2)') this%mass2
245 WRITE (excent_str,
'(ES10.2)') this%excent
246 WRITE (sema_str,
'(ES10.2)') this%semaaxis
247 CALL this%Info(
" primary mass: " // trim(mass1_str))
248 CALL this%Info(
" secondary mass: " // trim(mass2_str))
249 CALL this%Info(
" excentricity: " // trim(excent_str))
250 CALL this%Info(
" semi major axis: " // trim(sema_str))
251 IF(this%period.NE.0.0)
THEN
252 WRITE(omega_str,
'(ES10.2)') 2.0*pi/this%period
253 CALL this%Info(
" pattern speed: " // trim(omega_str))
256 IF (this%excent.LT.0.OR.this%excent.GE.1) &
257 CALL this%Error(
"InitGravity_binary",
"excentricity out of range, should be 0 <= e < 1 => check setup")
258 IF (this%semaaxis.LE.0) &
259 CALL this%Error(
"InitGravity_binary",
"semi major axis out of range, should be 0 < a => check setup")
269 TYPE(
dict_typ),
POINTER :: config,IO
274 CALL this%gravity_pointmass%SetOutput(mesh,physics,config,io)
277 CALL getattr(config,
"output/binpos", valwrite, 0)
278 IF (valwrite .EQ. 1) &
279 CALL setattr(io,
"binpos", this%pos)
301 REAL,
INTENT(IN) :: time,dt
307 CALL this%UpdatePositions(time)
311 this%posvec_prim_tmp(:,:,:,1) = this%pos(1,1)
312 this%posvec_prim_tmp(:,:,:,2) = this%pos(2,1)
313 this%posvec_prim_tmp(:,:,:,3) = this%pos(3,1)
314 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,
east,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,1,:))
315 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,
north,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,2,:))
316 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,
top,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,3,:))
317 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_prim_tmp,this%posvec_prim)
319 this%posvec_sec_tmp(:,:,:,1) = this%pos(1,2)
320 this%posvec_sec_tmp(:,:,:,2) = this%pos(2,2)
321 this%posvec_sec_tmp(:,:,:,3) = this%pos(3,2)
322 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,
east,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,1,:))
323 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,
north,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,2,:))
324 CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,
top,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,3,:))
325 CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_sec_tmp,this%posvec_sec)
331 this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_prim(:,:,:,:)
332 this%fposvec_prim(:,:,:,1,:) = mesh%posvec%faces(:,:,:,
east,:) - this%fposvec_prim(:,:,:,1,:)
333 this%fposvec_prim(:,:,:,2,:) = mesh%posvec%faces(:,:,:,
north,:) - this%fposvec_prim(:,:,:,2,:)
334 this%fposvec_prim(:,:,:,3,:) = mesh%posvec%faces(:,:,:,
top,:) - this%fposvec_prim(:,:,:,3,:)
336 this%posvec_sec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_sec(:,:,:,:)
337 this%fposvec_sec(:,:,:,1,:) = mesh%posvec%faces(:,:,:,
east,:) - this%fposvec_sec(:,:,:,1,:)
338 this%fposvec_sec(:,:,:,2,:) = mesh%posvec%faces(:,:,:,
north,:) - this%fposvec_sec(:,:,:,2,:)
339 this%fposvec_sec(:,:,:,3,:) = mesh%posvec%faces(:,:,:,
top,:) - this%fposvec_sec(:,:,:,3,:)
343 this%r_prim(:,:,:) = sqrt(this%posvec_prim(:,:,:,1)**2 &
344 + this%posvec_prim(:,:,:,2)**2 + this%posvec_prim(:,:,:,3)**2)
345 this%fr_prim(:,:,:,1) = sqrt(this%fposvec_prim(:,:,:,1,1)**2 &
346 + this%fposvec_prim(:,:,:,1,2)**2 + this%fposvec_prim(:,:,:,1,3)**2)
347 this%fr_prim(:,:,:,2) = sqrt(this%fposvec_prim(:,:,:,2,1)**2 &
348 + this%fposvec_prim(:,:,:,2,2)**2 + this%fposvec_prim(:,:,:,2,3)**2)
349 this%fr_prim(:,:,:,3) = sqrt(this%fposvec_prim(:,:,:,3,1)**2 &
350 + this%fposvec_prim(:,:,:,3,2)**2 + this%fposvec_prim(:,:,:,3,3)**2)
351 this%r_sec(:,:,:) = sqrt(this%posvec_sec(:,:,:,1)**2 &
352 + this%posvec_sec(:,:,:,2)**2 + this%posvec_sec(:,:,:,3)**2)
353 this%fr_sec(:,:,:,1) = sqrt(this%fposvec_sec(:,:,:,1,1)**2 &
354 + this%fposvec_sec(:,:,:,1,2)**2 + this%fposvec_sec(:,:,:,1,3)**2)
355 this%fr_sec(:,:,:,2) = sqrt(this%fposvec_sec(:,:,:,2,1)**2 &
356 + this%fposvec_sec(:,:,:,2,2)**2 + this%fposvec_sec(:,:,:,2,3)**2)
357 this%fr_sec(:,:,:,3) = sqrt(this%fposvec_sec(:,:,:,3,1)**2 &
358 + this%fposvec_sec(:,:,:,3,2)**2 + this%fposvec_sec(:,:,:,3,3)**2)
362 gm1 = physics%Constants%GN*this%GetMass_primary(time)
363 this%omega2%data4d(:,:,:,1) = gm1 / (this%r_prim(:,:,:)**3 + this%eps1**3)
366 gm2 = physics%Constants%GN*this%GetMass_secondary(time)
367 this%omega2%data4d(:,:,:,2) = gm2 / (this%r_sec(:,:,:)**3 + this%eps2**3)
370 CALL this%gravity_pointmass%CalcPotential(mesh,physics,this%mass,this%r_prim,this%fr_prim,this%pot_prim%data4d)
371 CALL this%gravity_pointmass%CalcPotential(mesh,physics,this%mass2,this%r_sec,this%fr_sec,this%pot_sec%data4d)
373 this%pot%data1d(:) = this%pot_prim%data1d(:) + this%pot_sec%data1d(:)
377 DO k=mesh%KGMIN,mesh%KGMAX
378 DO j=mesh%JGMIN,mesh%JGMAX
380 DO i=mesh%IGMIN,mesh%IGMAX
381 this%accel%data4d(i,j,k,1:physics%VDIM) = -this%omega2%data4d(i,j,k,1) * this%posvec_prim(i,j,k,1:physics%VDIM)&
382 -this%omega2%data4d(i,j,k,2) * this%posvec_sec(i,j,k,1:physics%VDIM)
404 REAL,
INTENT(IN) :: time
409 CALL kepler(time,this%period,this%excent,this%semaaxis,r,phi)
410 IF(r.LT.0)
CALL this%Error(
"UpdatePositions_binary",
"Solving Kepler-Equation failed!")
413 phi = phi - this%omega_rot * time
418 r1 = this%mass/(this%mass+this%mass2)*r
419 this%pos(1,2) = r1*cos(phi)
420 this%pos(2,2) = r1*sin(phi)
423 this%pos(:,1) = -this%pos(:,2) * this%mass2/this%mass
425 this%pos(:,1) = this%pos(:,1) + this%r0(:)
426 this%pos(:,2) = this%pos(:,2) + this%r0(:)
447 TYPE(
marray_base),
INTENT(INOUT) :: bccsound,h_ext,height
450 this%omega%data1d(:) = sqrt(
this%omega2%data2d(:,1) +
this%omega2%data2d(:,2))
461 REAL,
INTENT(IN) :: time
465 IF(time.LE.
this%switchon)
THEN
466 mass =
this%mass * sin(0.5*pi*time/
this%switchon)**2
477 REAL,
INTENT(IN) :: time
481 IF(time.LE.
this%switchon2)
THEN
482 mass =
this%mass2 * sin(0.5*pi*time/
this%switchon2)**2
494 DEALLOCATE(this%r0,this%r_sec,this%posvec_sec,this%posvec_sec_tmp,&
495 this%pot_prim,this%pot_sec,this%fr_sec,this%fposvec_sec,this%mass2)
497 CALL this%Finalize_base()
556 PURE SUBROUTINE kepler(time,period,excent,semax,r,phi)
560 REAL,
INTENT(IN) :: time
561 REAL,
INTENT(IN) :: period
562 REAL,
INTENT(IN) :: excent
563 REAL,
INTENT(IN) :: semax
564 REAL,
INTENT(OUT):: r
565 REAL,
INTENT(OUT):: phi
567 REAL,
DIMENSION(2):: plist
572 tau = 2.*pi*modulo(time,period)/period
575 IF ((tau.GE.0.0) .AND. (tau.LE.pi))
THEN
577 CALL getroot(
func,max(0.0,tau),min(pi,excent+tau),e,err,plist)
580 CALL getroot(
func,max(pi,tau-excent),min(2.*pi,tau),e,err,plist)
591 r = semax * (1.0 - excent*cose)
596 IF (cose.NE.-1.0)
THEN
597 phi = 2.0*atan(sign(sqrt(((1.+excent)*(1.-cose)) &
598 /((1.-excent)*(1.+cose))),pi-e))
606 PURE SUBROUTINE func(x,fx,plist)
609 REAL,
INTENT(IN) :: x
610 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
612 REAL,
INTENT(OUT) :: fx
614 fx = x - plist(1)*sin(x) - plist(2)
Dictionary for generic data types.
type(logging_base), save this
base module for numerical flux functions
generic gravity terms module providing functionaly common to all gravity terms
source terms module for gravitational acceleration due to two pointmasses
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
calculates the resultant gravitational accerleration of a binary system
real function getmass_primary(this, time)
get mass of the primary component, eventually scaled by switch on factor
subroutine infogravity(this)
write binary parameters to screen
real function getmass_secondary(this, time)
get mass of the primary component, eventually scaled by switch on factor
subroutine updatepositions(this, time)
compute new positions for both components of the binary system
pure subroutine func(x, fx, plist)
find root of the Kepler equation
pure subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
compute the pressure scale height of a disk for the binary potential
subroutine initgravity_binary(this, Mesh, Physics, config, IO)
Constructor of gravity binary module.
pure subroutine, public kepler(time, period, excent, semax, r, phi)
solve Kepler-Equation and return new distance r and position angle phi
source terms module for gravitational acceleration due to a point mass at the center of the coordinat...
elemental real function, public getdiskheight(cs, omega)
return pressure scale height of a geometrically thin Keplerian disk
subroutine finalize(this)
Destructor of logging_base class.
base class for mesh arrays
derived class for compound of mesh arrays
integer, parameter east
named constant for eastern boundary
subroutine setoutput(this, config, IO)
Setup mesh fields for i/o.
integer, parameter top
named constant for top boundary
integer, parameter north
named constant for northern boundary
physics module for 1D,2D and 3D non-isothermal Euler equations
physics module for 1D,2D and 3D isothermal Euler equations