49 PURE SUBROUTINE funcd(x,fx,dfx,plist)
52 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
53 REAL,
INTENT(OUT) :: fx,dfx
58 REAL,
PARAMETER :: msun = 1.989e+30
59 REAL,
PARAMETER :: gn = 6.67408e-11
61 REAL,
PARAMETER :: tsim = 20.0
62 REAL,
PARAMETER :: accmass = 1.0*msun
63 REAL,
PARAMETER :: gamma = 1.4
65 REAL,
PARAMETER :: rhoinf = 1.0e-20
66 REAL,
PARAMETER :: csinf = 1.0e+04
68 INTEGER,
PARAMETER :: mgeo = cylindrical
69 INTEGER,
PARAMETER :: xres = 100
70 INTEGER,
PARAMETER :: yres = 6
71 INTEGER,
PARAMETER :: zres = 1
72 REAL,
PARAMETER :: rin = 0.1
73 REAL,
PARAMETER :: rout = 2.0
74 REAL,
PARAMETER :: gpar = 1.0
76 INTEGER,
PARAMETER :: onum = 10
77 CHARACTER(LEN=256),
PARAMETER &
79 CHARACTER(LEN=256),
PARAMETER &
85 CLASS(
fosite),
ALLOCATABLE :: sim
86 REAL,
DIMENSION(:),
ALLOCATABLE :: sigma
87 REAL :: sum_numer, sum_denom
88 INTEGER :: n,den,vel,pre
96 IF (sim%GetRank().EQ.0)
THEN
107 CALL initdata(sim%Mesh, sim%Physics, sim%Fluxes, sim%Timedisc)
110 ok = .NOT.sim%aborted
112 IF (
ASSOCIATED(sim%Timedisc%solution))
THEN
113 ALLOCATE(sigma(sim%Physics%VNUM))
114 DO n=1,sim%Physics%VNUM
117 sum_numer = sum(abs(sim%Timedisc%pvar%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
118 sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n) &
119 -sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
120 sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
121 sum_denom = sum(abs(sim%Timedisc%solution%data4d(sim%Mesh%IMIN:sim%Mesh%IMAX,&
122 sim%Mesh%JMIN:sim%Mesh%JMAX,sim%Mesh%KMIN:sim%Mesh%KMAX,n)))
124 IF (sim%GetRank().GT.0)
THEN
125 CALL mpi_reduce(sum_numer,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
126 CALL mpi_reduce(sum_denom,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
128 CALL mpi_reduce(mpi_in_place,sum_numer,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
129 CALL mpi_reduce(mpi_in_place,sum_denom,1,default_mpi_real,mpi_sum,0,mpi_comm_world,sim%ierror)
131 sigma(n) = sum_numer / sum_denom
140 den = sim%Physics%DENSITY
141 vel = sim%Physics%XVELOCITY
142 pre = sim%Physics%PRESSURE
145 IF (sim%GetRank().EQ.0)
THEN
147tap_check(ok,
"stoptime reached")
149tap_check_small(sigma(den),1.0e-02,
"density deviation < 1%")
150tap_check_small(sigma(vel),1.0e-02,
"radial velocity deviation < 1%")
152tap_check_small(sigma(pre),1.0e-02,
"pressure deviation < 1%")
158 IF (
ALLOCATED(sigma))
DEALLOCATE(sigma)
167 CLASS(
fosite),
INTENT(INOUT) :: Sim
168 TYPE(dict_typ),
POINTER :: config
172 TYPE(dict_typ),
POINTER :: mesh, physics, boundary, datafile, &
173 grav, pmass, timedisc, fluxes
174 REAL :: x1,x2,y1,y2,z1,z2
177 rb = gn * accmass / csinf**2
189 bc(west) = no_gradients
193 bc(bottom) = no_gradients
194 bc(top) = no_gradients
196 CALL sim%Physics%Error(
"bondi2d::MakeConfig",
"geometry currently not supported")
200 "meshtype" / midpoint, &
211 "gparam" / (gpar*rb))
215 "western" / bc(west), &
216 "eastern" / bc(east), &
217 "southern" / bc(south), &
218 "northern" / bc(north), &
219 "bottomer" / bc(bottom), &
231 "variables" / conservative, &
232 "limiter" / monocent, &
237 "gtype" / pointmass, &
248 "method" / modified_euler, &
251 "stoptime" / (tsim * tau), &
252 "dtlimit" / (1.0e-6 * tau), &
253 "output/solution" / 1, &
258 "fileformat" / vtk, &
259 "filename" / (trim(odir) // trim(ofname)), &
264 "physics" / physics, &
265 "boundary" / boundary, &
267 "sources/grav" / grav, &
268 "timedisc" / timedisc, &
269 "datafile" / datafile)
276 CLASS(physics_base),
INTENT(IN) :: Physics
277 CLASS(mesh_base),
INTENT(IN) :: Mesh
278 CLASS(fluxes_base),
INTENT(IN) :: Fluxes
279 CLASS(timedisc_base),
INTENT(INOUT) :: Timedisc
284 CHARACTER(LEN=64) :: info_str
287 SELECT TYPE(pvar => timedisc%pvar)
288 TYPE IS(statevector_euler)
290 pvar%density%data1d(:) = rhoinf
291 pvar%velocity%data1d(:) = 0.
292 pvar%pressure%data1d(:) = rhoinf * csinf**2 / gamma
294 CALL physics%Error(
"bondi2d::InitData",
"only non-isothermal HD supported")
297 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
301 IF (
ASSOCIATED(timedisc%solution))
THEN
303 SELECT TYPE(pvar => timedisc%solution)
304 TYPE IS(statevector_euler)
305 DO k=mesh%KMIN,mesh%KMAX
306 DO j=mesh%JMIN,mesh%JMAX
307 DO i=mesh%IMIN,mesh%IMAX
308 r = mesh%radius%bcenter(i,j,k)
309 CALL bondi(r/rb,gamma,rhoinf,csinf,rho,vr)
310 cs2 = csinf**2 * (rho/rhoinf)**(gamma-1.0)
311 pvar%density%data3d(i,j,k) = rho
316 pvar%velocity%data4d(i,j,k,l) = sign(&
317 max(epsilon(cs2),abs(vr*mesh%posvec%bcenter(i,j,k,l)/r)), &
318 vr*mesh%posvec%bcenter(i,j,k,l))
320 pvar%pressure%data3d(i,j,k) = rho * cs2 / gamma
329 SELECT TYPE(beast => timedisc%Boundary%boundary(east)%p)
330 CLASS IS (boundary_custom)
332 CALL beast%SetCustomBoundaries(mesh,physics, &
333 (/custom_fixed,custom_nograd,custom_fixed,custom_fixed/))
334 DO k=mesh%KMIN,mesh%KMAX
335 DO j=mesh%JMIN,mesh%JMAX
338 r = mesh%radius%bcenter(mesh%IMAX+i,j,k)
339 CALL bondi(r/rb,gamma,rhoinf,csinf,rho,vr)
340 cs2 = csinf**2 * (rho/rhoinf)**(gamma-1.0)
341 beast%data(i,j,k,physics%DENSITY) = rho
343 beast%data(i,j,k,physics%XVELOCITY+l-1) &
344 = vr*mesh%posvec%bcenter(mesh%IMAX+i,j,k,l)/r
346 beast%data(i,j,k,physics%PRESSURE) = rho * cs2 / gamma
353 CALL mesh%Info(
" DATA-----> initial condition: 2D Bondi accretion")
354 WRITE(info_str,
"(ES11.3)") rb
355 CALL mesh%Info(
" " //
"Bondi radius: " &
356 // trim(info_str) //
" m")
357 WRITE(info_str,
"(ES11.3)") tau
358 CALL mesh%Info(
" " //
"Free fall time: " &
359 // trim(info_str) //
" s")
362 SUBROUTINE bondi(r,gamma,rhoinf,csinf,rho,vr)
378 REAL,
INTENT(IN) :: r,gamma,rhoinf,csinf
379 REAL,
INTENT(OUT) :: rho,vr
381 REAL,
PARAMETER :: xacc = 1.0e-6
382 REAL :: gp1,gm1,rc,chi,lambda,psi,gr
391 rc = (3.0-gamma) / 2.0
394 lambda = rc**(-rc/gm1)
398 gr = chi**(2.*gm1/gp1) * (1./gm1 + 1./r)
408 rho = rhoinf * chi**(-2./gp1) / psi
409 vr = -csinf * chi**(-gm1/gp1) * psi
415PURE SUBROUTINE funcd(y,fy,dfy,plist)
418 REAL,
INTENT(IN) :: y
419 REAL,
INTENT(IN),
DIMENSION(:),
OPTIONAL :: plist
420 REAL,
INTENT(OUT) :: fy,dfy
425 fy = 0.5*y*y + y**(-plist(1)) / plist(1) - plist(2)
426 dfy = y - y**(-plist(1)-1.)
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
pure subroutine funcd(y, fy, dfy, plist)
subroutine makeconfig(Sim, config)
program bondi2d
Program and data initialization for 2D Bondi accretion.
subroutine bondi(r, gamma, rhoinf, csinf, rho, vr)
subroutine, public getroot_newton(funcd, x1, x2, root, error, plist, xm, iterations)