gravity_binary.f90
Go to the documentation of this file.
1 !#############################################################################
2 !# #
3 !# fosite - 3D hydrodynamical simulation program #
4 !# module: gravity_binary.f90 #
5 !# #
6 !# Copyright (C) 2010-2019 #
7 !# Tobias Illenseer <tillense@astrophysik.uni-kiel.de> #
8 !# Björn Sperling <sperling@astrophysik.uni-kiel.de> #
9 !# Anna Feiler <afeiler@astrophysik.uni-kiel.de> #
10 !# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
11 !# #
12 !# This program is free software; you can redistribute it and/or modify #
13 !# it under the terms of the GNU General Public License as published by #
14 !# the Free Software Foundation; either version 2 of the License, or (at #
15 !# your option) any later version. #
16 !# #
17 !# This program is distributed in the hope that it will be useful, but #
18 !# WITHOUT ANY WARRANTY; without even the implied warranty of #
19 !# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
20 !# NON INFRINGEMENT. See the GNU General Public License for more #
21 !# details. #
22 !# #
23 !# You should have received a copy of the GNU General Public License #
24 !# along with this program; if not, write to the Free Software #
25 !# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
26 !# #
27 !#############################################################################
40 !----------------------------------------------------------------------------!
52 !----------------------------------------------------------------------------!
53 MODULE gravity_binary_mod
54  USE gravity_pointmass_mod
55  USE gravity_base_mod
56  USE fluxes_base_mod
57  USE physics_base_mod
58  USE mesh_base_mod
59  USE logging_base_mod
60  USE marray_base_mod
61  USE marray_compound_mod
62  USE common_dict
63 #ifdef PARALLEL
64 #ifdef HAVE_MPI_MOD
65  USE mpi
66 #endif
67 #endif
68  IMPLICIT NONE
69 #ifdef PARALLEL
70 #ifdef HAVE_MPIF_H
71  include 'mpif.h'
72 #endif
73 #endif
74  !--------------------------------------------------------------------------!
75  PRIVATE
76  TYPE, EXTENDS(gravity_pointmass) :: gravity_binary
77  REAL, DIMENSION(:,:,:), POINTER :: r_sec
78  REAL, DIMENSION(:,:,:,:), POINTER :: posvec_sec
79  REAL, DIMENSION(:,:,:,:), POINTER :: posvec_sec_tmp
80  REAL, DIMENSION(:,:,:,:), POINTER :: fr_sec
81  REAL, DIMENSION(:,:,:,:,:), POINTER :: fposvec_sec
82  REAL, POINTER :: mass2
83  REAL :: excent
84  REAL :: semaaxis
85  REAL :: period
86  REAL :: eps1,eps2
87  REAL :: switchon2
88  REAL :: omega_rot
89  REAL, DIMENSION(:,:,:,:), POINTER :: pot_sec
90  CONTAINS
91  PROCEDURE :: initgravity_binary
92  PROCEDURE :: updategravity_single
93  PROCEDURE :: updatepositions
94  PROCEDURE :: infogravity
95  PROCEDURE :: setoutput
96  PROCEDURE :: getmass_primary
97  PROCEDURE :: getmass_secondary
98  PROCEDURE :: calcdiskheight_single
99  PROCEDURE :: finalize
100  END TYPE
101 
102  PUBLIC :: &
103  ! types
104  gravity_binary
105  !--------------------------------------------------------------------------!
106 
107 CONTAINS
113  SUBROUTINE initgravity_binary(this,Mesh,Physics,config,IO)
114  USE physics_euler_mod, ONLY: physics_euler
115  USE physics_eulerisotherm_mod, ONLY: physics_eulerisotherm
116  IMPLICIT NONE
117  !------------------------------------------------------------------------!
118  CLASS(gravity_binary), INTENT(INOUT) :: this
119  CLASS(mesh_base), INTENT(IN) :: Mesh
120  CLASS(physics_base), INTENT(IN) :: Physics
121  TYPE(Dict_TYP),POINTER :: config
123  TYPE(Dict_TYP),POINTER :: IO
124  !------------------------------------------------------------------------!
125  INTEGER :: mesh_rot
126  INTEGER :: err
127  !------------------------------------------------------------------------!
128  ! init base class
129  CALL this%InitGravity(mesh,physics,"binary point masses",config,io)
130 
131  SELECT TYPE(physics)
132  TYPE IS(physics_euler)
133  ! do nothing
134  TYPE IS(physics_eulerisotherm)
135  ! do nothing
136  CLASS DEFAULT
137  CALL this%Error("InitGravity_binary", &
138  "Physics not supported for binary gravity term")
139  END SELECT
140 
141  ALLOCATE(&
142  this%mass,this%mass2,this%pos(3,2),this%r0(3), &
143  this%omega2(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,2), &
144  this%omega(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX), &
145  this%r_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
146  this%r_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX),&
147  this%posvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
148  this%posvec_prim_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
149  this%posvec_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
150  this%posvec_sec_tmp(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3),&
151  this%pot(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,4), &
152  this%pot_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,4), &
153  this%pot_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,4), &
154  this%fr_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
155  this%fr_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3), &
156  this%fposvec_prim(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),&
157  this%fposvec_sec(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX,3,3),& ! last two entries (EAST,NORTH)x(dim1,dim2)
158  stat = err)
159  IF (err.NE.0) CALL this%Error("InitGravity_pointmass", "Unable allocate memory!")
160 
161  ! mass of primary component
162  CALL getattr(config, "mass1", this%mass, 1.0)
163  CALL setattr(io, "mass1", ref(this%mass))
164 
165  ! mass of secondary component
166  CALL getattr(config, "mass2", this%mass2, 1.0)
167  CALL setattr(io, "mass2", ref(this%mass2))
168 
169  ! excentricity
170  CALL getattr(config, "excentricity", this%excent, 0.0)
171 
172  ! semi major axis
173  CALL getattr(config, "semimayoraxis", this%semaaxis, 1.0)
174 
175  ! Softening parameter
176  CALL getattr(config, "softening1", this%eps1, 0.0)
177 
178  CALL getattr(config, "softening2", this%eps2, 0.0)
179 
180  ! soft switch on for primary component
181  CALL getattr(config, "switchon1", this%switchon, -1.0)
182 
183  ! soft switch on for secondary component
184  CALL getattr(config, "switchon2", this%switchon2, -1.0)
185 
186  ! rotating reference frame
187  CALL getattr(config, "omega_rot", this%omega_rot, 0.0)
188 
189  ! check whether rotating mesh is activated
190  CALL getattr(config, "mesh_rot", mesh_rot, 0)
191  IF (mesh_rot.EQ.1) THEN
192  this%omega_rot = mesh%OMEGA
193  this%r0(:) = mesh%rotcent(:)
194  IF (mesh%omega.EQ.0.0) &
195  CALL this%Warning("InitGravity_binary","Rotating mesh enabled but OMEGA set to zero in mesh")
196  ELSE
197  this%r0(:) = 0.0
198  END IF
199 
200  ! reset mass flux
201  this%mdot = 0.0
202 
203  ! set period
204  this%period = 2.*pi*sqrt(this%semaaxis/(this%mass+this%mass2)/physics%constants%GN)*this%semaaxis
205 
206  ! reset omega and omega**2
207  this%omega = 0.0
208  this%omega2 = 0.0
209 
210  ! reset potentials
211  this%pot = 0.0
212  this%pot_prim = 0.0
213  this%pot_sec = 0.0
214 
215  ! set ghost cell data to one to avoid zero division
216  this%omega(mesh%IGMIN:mesh%IMIN+mesh%IM1,:,:) = 1.0
217  this%omega(mesh%IMAX+mesh%IP1:mesh%IGMAX,:,:) = 1.0
218  this%omega(:,mesh%JGMIN:mesh%JMIN+mesh%JM1,:) = 1.0
219  this%omega(:,mesh%JMAX+mesh%JP1:mesh%JGMAX,:) = 1.0
220  this%omega(:,:,mesh%KGMIN:mesh%KMIN+mesh%KM1) = 1.0
221  this%omega(:,:,mesh%KMAX+mesh%KP1:mesh%KGMAX) = 1.0
222 
223  !initialise output
224  CALL this%SetOutput(mesh,physics,config,io)
225  CALL this%InfoGravity()
226 
227  END SUBROUTINE initgravity_binary
228 
229 
231  SUBROUTINE infogravity(this)
232  IMPLICIT NONE
233  !------------------------------------------------------------------------!
234  CLASS(gravity_binary), INTENT(IN) :: this
235  !------------------------------------------------------------------------!
236  CHARACTER(LEN=32) :: mass1_str,mass2_str,excent_str,sema_str,omega_str
237  !------------------------------------------------------------------------!
238  WRITE (mass1_str,'(ES8.2)') this%mass
239  WRITE (mass2_str,'(ES8.2)') this%mass2
240  WRITE (excent_str,'(ES8.2)') this%excent
241  WRITE (sema_str,'(ES8.2)') this%semaaxis
242  CALL this%Info(" primary mass: " // trim(mass1_str))
243  CALL this%Info(" secondary mass: " // trim(mass2_str))
244  CALL this%Info(" excentricity: " // trim(excent_str))
245  CALL this%Info(" semi major axis: " // trim(sema_str))
246  IF(this%period.NE.0.0) THEN
247  WRITE(omega_str,'(ES8.2)') 2.0*pi/this%period
248  CALL this%Info(" pattern speed: " // trim(omega_str))
249  END IF
250  END SUBROUTINE infogravity
251 
252 
253  SUBROUTINE setoutput(this,Mesh,Physics,config,IO)
254  IMPLICIT NONE
255  !------------------------------------------------------------------------!
256  CLASS(gravity_binary), INTENT(INOUT) :: this
257  CLASS(mesh_base), INTENT(IN) :: Mesh
258  CLASS(physics_base), INTENT(IN) :: Physics
259  TYPE(Dict_TYP), POINTER :: config,IO
260  !------------------------------------------------------------------------!
261  INTEGER :: valwrite
262  !------------------------------------------------------------------------!
263 
264  CALL this%gravity_pointmass%SetOutput(mesh,physics,config,io)
265 
266  ! output cartesian positions of both components of the binary system
267  CALL getattr(config, "output/binpos", valwrite, 0)
268  IF (valwrite .EQ. 1) &
269  CALL setattr(io, "binpos", this%pos)
270 
271  END SUBROUTINE setoutput
272 
273 
283  SUBROUTINE updategravity_single(this,Mesh,Physics,Fluxes,pvar,time,dt)
284  IMPLICIT NONE
285  !------------------------------------------------------------------------!
286  CLASS(gravity_binary), INTENT(INOUT) :: this
287  CLASS(mesh_base), INTENT(IN) :: Mesh
288  CLASS(physics_base), INTENT(IN) :: Physics
289  CLASS(fluxes_base), INTENT(IN) :: Fluxes
290  CLASS(marray_compound),INTENT(INOUT) :: pvar
291  REAL, INTENT(IN) :: time,dt
292  !------------------------------------------------------------------------!
293  REAL :: GM1, GM2
294  INTEGER :: i,j,k
295  !------------------------------------------------------------------------!
296  ! update the positions of the binary system
297  CALL this%UpdatePositions(time)
298 
299  ! compute curvilinear components of the position vectors
300  ! 1. primary component
301  this%posvec_prim_tmp(:,:,:,1) = this%pos(1,1)
302  this%posvec_prim_tmp(:,:,:,2) = this%pos(2,1)
303  this%posvec_prim_tmp(:,:,:,3) = this%pos(3,1)
304  CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,east,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,1,:))
305  CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,north,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,2,:))
306  CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,top,:),this%posvec_prim_tmp(:,:,:,:),this%fposvec_prim(:,:,:,3,:))
307  CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_prim_tmp,this%posvec_prim)
308  ! 2. secondary component
309  this%posvec_sec_tmp(:,:,:,1) = this%pos(1,2)
310  this%posvec_sec_tmp(:,:,:,2) = this%pos(2,2)
311  this%posvec_sec_tmp(:,:,:,3) = this%pos(3,2)
312  CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,east,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,1,:))
313  CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,north,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,2,:))
314  CALL mesh%Geometry%Convert2Curvilinear(mesh%curv%faces(:,:,:,top,:),this%posvec_sec_tmp(:,:,:,:),this%fposvec_sec(:,:,:,3,:))
315  CALL mesh%Geometry%Convert2Curvilinear(mesh%bcenter,this%posvec_sec_tmp,this%posvec_sec)
316 
317  ! subtract the result from the position vectors of all cell bary centers
318  ! this gives you the curvilinear components of all vectors pointing
319  ! from the point mass of each component of the binary system to the bary
320  ! centers of all cells on the mesh
321  this%posvec_prim(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_prim(:,:,:,:)
322  this%fposvec_prim(:,:,:,1,:) = mesh%posvec%faces(:,:,:,east,:) - this%fposvec_prim(:,:,:,1,:) ! EAST
323  this%fposvec_prim(:,:,:,2,:) = mesh%posvec%faces(:,:,:,north,:) - this%fposvec_prim(:,:,:,2,:) ! NORTH
324  this%fposvec_prim(:,:,:,3,:) = mesh%posvec%faces(:,:,:,top,:) - this%fposvec_prim(:,:,:,3,:) ! TOP
325 
326  this%posvec_sec(:,:,:,:) = mesh%posvec%bcenter(:,:,:,:) - this%posvec_sec(:,:,:,:)
327  this%fposvec_sec(:,:,:,1,:) = mesh%posvec%faces(:,:,:,east,:) - this%fposvec_sec(:,:,:,1,:) ! EAST
328  this%fposvec_sec(:,:,:,2,:) = mesh%posvec%faces(:,:,:,north,:) - this%fposvec_sec(:,:,:,2,:) ! NORTH
329  this%fposvec_sec(:,:,:,3,:) = mesh%posvec%faces(:,:,:,top,:) - this%fposvec_sec(:,:,:,3,:) ! TOP
330 
331  ! compute the distances between each component of the binary system
332  ! and all cell bary centers
333  this%r_prim(:,:,:) = sqrt(this%posvec_prim(:,:,:,1)**2+this%posvec_prim(:,:,:,2)**2+this%posvec_prim(:,:,:,3)**2)
334  this%fr_prim(:,:,:,1) = sqrt(this%fposvec_prim(:,:,:,1,1)**2+this%fposvec_prim(:,:,:,1,2)**2+this%fposvec_prim(:,:,:,1,3)**2) ! shifted EAST-faces
335  this%fr_prim(:,:,:,2) = sqrt(this%fposvec_prim(:,:,:,2,1)**2+this%fposvec_prim(:,:,:,2,2)**2+this%fposvec_prim(:,:,:,2,3)**2) ! shifted NORTH-faces
336  this%fr_prim(:,:,:,3) = sqrt(this%fposvec_prim(:,:,:,3,1)**2+this%fposvec_prim(:,:,:,3,2)**2+this%fposvec_prim(:,:,:,3,3)**2) ! shifted TOP-faces
337  this%r_sec(:,:,:) = sqrt(this%posvec_sec(:,:,:,1)**2+this%posvec_sec(:,:,:,2)**2+this%posvec_sec(:,:,:,3)**2)
338  this%fr_sec(:,:,:,1) = sqrt(this%fposvec_sec(:,:,:,1,1)**2+this%fposvec_sec(:,:,:,1,2)**2+this%fposvec_sec(:,:,:,1,3)**2) ! shifted EAST-faces
339  this%fr_sec(:,:,:,2) = sqrt(this%fposvec_sec(:,:,:,2,1)**2+this%fposvec_sec(:,:,:,2,2)**2+this%fposvec_sec(:,:,:,2,3)**2) ! shifted NORTH-faces
340  this%fr_sec(:,:,:,3) = sqrt(this%fposvec_sec(:,:,:,3,1)**2+this%fposvec_sec(:,:,:,3,2)**2+this%fposvec_sec(:,:,:,3,3)**2) ! shifted NORTH-faces
341 
342  ! compute square of Keplerian velocities for updated time value
343  ! with respect to PRIMARY star
344  gm1 = physics%Constants%GN*this%GetMass_primary(time)
345  this%omega2(:,:,:,1) = gm1 / (this%r_prim(:,:,:)**3 + this%eps1**3)
346 
347  ! and SECONDARY star
348  gm2 = physics%Constants%GN*this%GetMass_secondary(time)
349  this%omega2(:,:,:,2) = gm2 / (this%r_sec(:,:,:)**3 + this%eps2**3)
350 
351  ! potential calculation for binary and use pointmassroutine for it
352  CALL this%gravity_pointmass%CalcPotential(mesh,physics,this%mass,this%r_prim,this%fr_prim,this%pot_prim)
353  CALL this%gravity_pointmass%CalcPotential(mesh,physics,this%mass2,this%r_sec,this%fr_sec,this%pot_sec)
354 
355  this%pot(:,:,:,:) = this%pot_prim(:,:,:,:) + this%pot_sec(:,:,:,:)
356 
357  ! set curvilinear components of the gravitational acceleration
358 !NEC$ collapse
359  DO k=mesh%KGMIN,mesh%KGMAX
360  DO j=mesh%JGMIN,mesh%JGMAX
361 !NEC$ ivdep
362  DO i=mesh%IGMIN,mesh%IGMAX
363  this%accel%data4d(i,j,k,1:physics%VDIM) = -this%omega2(i,j,k,1) * this%posvec_prim(i,j,k,1:physics%VDIM)&
364  -this%omega2(i,j,k,2) * this%posvec_sec(i,j,k,1:physics%VDIM)
365  END DO
366  END DO
367  END DO
368 
369  END SUBROUTINE updategravity_single
370 
371 
436  SUBROUTINE updatepositions(this,time)
437  USE roots
438  IMPLICIT NONE
439  !------------------------------------------------------------------------!
440  CLASS(gravity_binary), INTENT(INOUT) :: this
441  REAL, INTENT(IN) :: time
442  !------------------------------------------------------------------------!
443  REAL :: E,cosE,r,phi,r1
444  REAL :: tau,excent
445  REAL, DIMENSION(2):: plist
446  INTEGER :: err
447  !------------------------------------------------------------------------!
448  ! mean anomaly
449  tau = 2.*pi*modulo(time,this%period)/this%period
450  excent = this%excent
451 
452  plist(1)=excent
453  plist(2)=tau
454  ! solve the Kepler equation: E-tau-excent*sin(E) = 0
455  ! i.e. compute the intersection of f(E) = E-tau and g(E) = excent*sin(E)
456  IF ((tau.GE.0.0) .AND. (tau.LE.pi)) THEN
457  ! intersection lies ABOVE abscissa
458  CALL getroot(func,max(0.0,tau),min(pi,excent+tau),e,err,plist)
459  ELSE
460  ! intersection lies BELOW abscissa
461  CALL getroot(func,max(pi,tau-excent),min(2.*pi,tau),e,err,plist)
462  END IF
463  IF(err.NE.0) &
464  CALL this%Error("UpdatePositions_binary","Error in root finding!")
465 
466  ! compute positions
467  cose = cos(e)
468  ! distance to the origin of primary star position
469  r = this%semaaxis * (1.0 - excent*cose)
470  ! position angle of primary star
471  ! phi1 = 2.0*ATAN(SQRT((1.0+excent)/(1.0-excent))*TAN(0.5*E))
472  ! A little bit more complicated but this avoids one evaluation of TAN;
473  ! uses TAN(x/2) = +/- SQRT((1-cos(x))/(1+cos(x)))
474  IF (cose.NE.-1.0) THEN
475  phi = 2.0*atan(sign(sqrt(((1.+excent)*(1.-cose)) &
476  /((1.-excent)*(1.+cose))),pi-e))
477  ELSE
478  ! since lim_{x->0} atan(q/x) = pi/2 for q>0
479  phi = pi
480  END IF
481 
482  !transform to rotating reference frame if necessary
483  phi = phi - this%omega_rot * time
484 
485  ! cartesian coordinates of SECONDARY component (right at t=0)
486  ! \todo This implementation only works in 2D
487  r1 = this%mass/(this%mass+this%mass2)*r
488  this%pos(1,2) = r1*cos(phi)
489  this%pos(2,2) = r1*sin(phi)
490  this%pos(3,2) = 0.0
491  ! and PRIMARY COMPONENT
492  this%pos(:,1) = -this%pos(:,2) * this%mass2/this%mass
493  ! shift with respect to center of rotation
494  this%pos(:,1) = this%pos(:,1) + this%r0(:)
495  this%pos(:,2) = this%pos(:,2) + this%r0(:)
496 
497  END SUBROUTINE updatepositions
498 
510  PURE SUBROUTINE calcdiskheight_single(this,Mesh,Physics,pvar,bccsound,h_ext,height)
511  IMPLICIT NONE
512  !------------------------------------------------------------------------!
513  CLASS(gravity_binary), INTENT(INOUT) :: this
514  CLASS(mesh_base), INTENT(IN) :: mesh
515  CLASS(physics_base), INTENT(IN) :: physics
516  CLASS(marray_compound), INTENT(INOUT) :: pvar
517  TYPE(marray_base), INTENT(INOUT) :: bccsound,h_ext,height
518  !------------------------------------------------------------------------!
519  ! compute the effective omega
520  this%omega(:,:,:) = sqrt(this%omega2(:,:,:,1) + this%omega2(:,:,:,2))
521  ! compute the disk height
522  height%data3d(:,:,:) = getdiskheight(bccsound%data3d(:,:,:),this%omega(:,:,:))
523  END SUBROUTINE calcdiskheight_single
524 
525 
527  FUNCTION getmass_primary(this,time) RESULT(mass)
528  IMPLICIT NONE
529  !------------------------------------------------------------------------!
530  CLASS(gravity_binary), INTENT(IN) :: this
531  REAL, INTENT(IN) :: time
532  !------------------------------------------------------------------------!
533  REAL :: mass
534  !------------------------------------------------------------------------!
535  IF(time.LE.this%switchon) THEN
536  mass = this%mass * sin(0.5*pi*time/this%switchon)**2
537  ELSE
538  mass = this%mass
539  END IF
540  END FUNCTION getmass_primary
541 
543  FUNCTION getmass_secondary(this,time) RESULT(mass)
544  IMPLICIT NONE
545  !------------------------------------------------------------------------!
546  CLASS(gravity_binary), INTENT(IN) :: this
547  REAL, INTENT(IN) :: time
548  !------------------------------------------------------------------------!
549  REAL :: mass
550  !------------------------------------------------------------------------!
551  IF(time.LE.this%switchon2) THEN
552  mass = this%mass2 * sin(0.5*pi*time/this%switchon2)**2
553  ELSE
554  mass = this%mass2
555  END IF
556  END FUNCTION getmass_secondary
557 
559  SUBROUTINE finalize(this)
560  IMPLICIT NONE
561  !------------------------------------------------------------------------!
562  CLASS(gravity_binary), INTENT(INOUT) :: this
563  !------------------------------------------------------------------------!
564  DEALLOCATE(this%mass,this%mass2,this%pos,this%r0,this%omega2,this%omega,&
565  this%r_prim,this%r_sec,this%posvec_prim,this%posvec_prim_tmp, &
566  this%posvec_sec,this%posvec_sec_tmp,this%pot,&
567  this%pot_prim,this%pot_sec,this%fr_prim,this%fr_sec,this%fposvec_prim, &
568  this%fposvec_sec)
569 
570  CALL this%Finalize_base()
571  END SUBROUTINE finalize
572 
573 
575  PURE SUBROUTINE func(x,fx,plist)
576  IMPLICIT NONE
577  !------------------------------------------------------------------------!
578  REAL, INTENT(IN) :: x
579  REAL, INTENT(IN), DIMENSION(:), OPTIONAL :: plist
580  !plist[1]=excent, plist[2]=tau
581  REAL, INTENT(OUT) :: fx
582  !------------------------------------------------------------------------!
583  fx = x - plist(1)*sin(x) - plist(2)
584  END SUBROUTINE func
585 
586 END MODULE gravity_binary_mod
mesh_base
logging_base_mod::finalize
subroutine finalize(this)
Destructor of common class.
Definition: logging_base.f90:155
gravity_binary_mod::infogravity
subroutine infogravity(this)
write binary parameters to screen
Definition: gravity_binary.f90:232
gravity_binary_mod::initgravity_binary
subroutine initgravity_binary(this, Mesh, Physics, config, IO)
Constructor of gravity binary module.
Definition: gravity_binary.f90:114
gravity_binary_mod::gravity_binary
Definition: gravity_binary.f90:76
common_dict::this
type(logging_base), save this
Definition: common_dict.f90:149
marray_compound_mod
derived class for compound of mesh arrays
Definition: marray_compound.f90:35
marray_base_mod
base class for mesh arrays
Definition: marray_base.f90:36
gravity_binary_mod::getmass_secondary
real function getmass_secondary(this, time)
get mass of the primary component, eventually scaled by switch on factor
Definition: gravity_binary.f90:544
physics_eulerisotherm_mod::physics_eulerisotherm
Definition: physics_eulerisotherm.f90:57
logging_base_mod
Basic fosite module.
Definition: logging_base.f90:36
gravity_binary_mod::updatepositions
subroutine updatepositions(this, time)
compute new positions for both components of the binary system
Definition: gravity_binary.f90:437
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_base_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: gravity_base.f90:47
gravity_base_mod::setoutput
subroutine setoutput(this, Mesh, Physics, config, IO)
Definition: gravity_base.f90:145
gravity_binary_mod::calcdiskheight_single
pure subroutine calcdiskheight_single(this, Mesh, Physics, pvar, bccsound, h_ext, height)
compute the pressure scale height of a disk for the binary potential
Definition: gravity_binary.f90:511
gravity_binary_mod::updategravity_single
subroutine updategravity_single(this, Mesh, Physics, Fluxes, pvar, time, dt)
calculates the resultant gravitational accerleration of a binary system
Definition: gravity_binary.f90:284
marray_base_mod::marray_base
basic mesh array class
Definition: marray_base.f90:64
physics_eulerisotherm_mod
physics module for 1D,2D and 3D isothermal Euler equations
Definition: physics_eulerisotherm.f90:44
physics_base_mod::physics_base
named integer constants for flavour of state vectors
Definition: physics_base.f90:69
physics_base_mod
Basic physics module.
Definition: physics_base.f90:55
marray_compound_mod::marray_compound
data types and methods
Definition: marray_compound.f90:46
roots
root finding subroutines
Definition: roots.f90:45
common_dict
Dictionary for generic data types.
Definition: common_dict.f90:61
gravity_pointmass_mod::gravity_pointmass
Definition: gravity_pointmass.f90:77
physics_euler_mod
physics module for 1D,2D and 3D non-isothermal Euler equations
Definition: physics_euler.f90:41
gravity_binary_mod::getmass_primary
real function getmass_primary(this, time)
get mass of the primary component, eventually scaled by switch on factor
Definition: gravity_binary.f90:528
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:465
physics_euler_mod::physics_euler
Definition: physics_euler.f90:54
fluxes_base_mod
base module for numerical flux functions
Definition: fluxes_base.f90:39
gravity_binary_mod::func
pure subroutine func(x, fx, plist)
find root of the Kepler equation
Definition: gravity_binary.f90:576
gravity_binary_mod
source terms module for gravitational acceleration due to two pointmasses
Definition: gravity_binary.f90:53
Generated by   doxygen 1.8.14