orbitingcylinders.f90
Go to the documentation of this file.
1!#############################################################################
2!# #
3!# fosite - 3D hydrodynamical simulation program #
4!# module: orbitingcylinders.f90 #
5!# #
6!# Copyright (C) 2013-2024 #
7!# Manuel Jung <mjung@astrophysik.uni-kiel.de> #
8!# Jannes Klee <jklee@astrophysik.uni-kiel.de> #
9!# #
10!# This program is free software; you can redistribute it and/or modify #
11!# it under the terms of the GNU General Public License as published by #
12!# the Free Software Foundation; either version 2 of the License, or (at #
13!# your option) any later version. #
14!# #
15!# This program is distributed in the hope that it will be useful, but #
16!# WITHOUT ANY WARRANTY; without even the implied warranty of #
17!# MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, GOOD TITLE or #
18!# NON INFRINGEMENT. See the GNU General Public License for more #
19!# details. #
20!# #
21!# You should have received a copy of the GNU General Public License #
22!# along with this program; if not, write to the Free Software #
23!# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. #
24!# #
25!#############################################################################
26
27!----------------------------------------------------------------------------!
152!----------------------------------------------------------------------------!
153PROGRAM orbitingcylinders
154 USE fosite_mod
155 IMPLICIT NONE
156 !--------------------------------------------------------------------------!
157 ! problem setup
158 REAL, PARAMETER :: sigma = 0.1
159 REAL, DIMENSION(2), PARAMETER :: x1 = (/ 1.0, 0. /)
160 REAL, DIMENSION(2), PARAMETER :: x2 = (/ 1.0, pi /)
161 REAL, DIMENSION(2), PARAMETER :: x3 = (/ 0.9, 3./4.*pi /)
162 REAL, DIMENSION(2), PARAMETER :: x4 = (/ 1.0, -pi/2. /)
163
164 ! general constants
165 REAL, PARAMETER :: gn = 1.
166 REAL, PARAMETER :: p0 = 1.
167 REAL, PARAMETER :: p = 10.
168
169 REAL, PARAMETER :: flaring = 0.05
170 REAL, PARAMETER :: rg = 8.31
171 REAL, PARAMETER :: mu = 6.02e-04
172 REAL, PARAMETER :: gamma = 5./3.
173
174 REAL, PARAMETER :: rmin = 0.4
175 REAL, PARAMETER :: rmax = 2.0
176 INTEGER, PARAMETER :: yres = 128
177 INTEGER, PARAMETER :: xres = yres/4
178 INTEGER, PARAMETER :: onum = p * 1
179 REAL, PARAMETER :: omega = 1.0
180 !--------------------------------------------------------------------------!
181 CLASS(fosite), ALLOCATABLE :: sim
182 INTEGER :: green
183 CHARACTER(LEN=1) :: str
184 REAL, DIMENSION(2) :: pot_err
185 REAL, DIMENSION(:,:,:), POINTER :: numpot => null(), anapot => null()
186 !--------------------------------------------------------------------------!
187
188 DO green=1,3
189 ALLOCATE(sim)
190 CALL sim%InitFosite()
191 ! basic configuration for all tests
192 CALL makeconfig(sim,sim%config)
193
194 ! set data file name and Green's function type
195 WRITE(str,'(I1)') green
196 CALL setattr(sim%config, "/datafile/filename",("orbitingcylinders-test" // str))
197 CALL setattr(sim%config, "/sources/grav/self/green",green)
198
199 SELECT CASE(green)
200 CASE(1,2)
201 CALL setattr(sim%config, "/timedisc/stoptime",1e-5*p0)
202 CALL setattr(sim%config, "/datafile/count",1)
203 CASE(3)
204 CALL setattr(sim%config, "/mesh/xmin", 0.2)
205 CALL setattr(sim%config, "/mesh/xmax", 1.8)
206 CALL setattr(sim%config, "/mesh/inum", 64)
207 CALL setattr(sim%config, "/mesh/jnum", 192)
208 CALL setattr(sim%config, "/timedisc/stoptime",p*p0)
209 CALL setattr(sim%config, "/datafile/count",onum)
210 END SELECT
211
212 CALL sim%Setup()
213 CALL initdata(sim%Timedisc,sim%Mesh,sim%Physics,sim%Fluxes,sim%Sources)
214 CALL sim%FirstStep()
215
216 SELECT CASE(green)
217 CASE(1,2)
218! CALL Sim%Datafile%WriteDataset(Sim%Mesh,Sim%Physics,Sim%Fluxes,&
219! Sim%Timedisc,Sim%config,Sim%IO)
220 ! compute max of local/global relative error
221 pot_err(1) = maxval(abs((numpot(:,:,:)-anapot(:,:,:))/anapot(:,:,:)), &
222 mask=sim%Mesh%without_ghost_zones%mask3d(:,:,:))
223 pot_err(2) = sum(abs(numpot-anapot),mask=sim%Mesh%without_ghost_zones%mask3d(:,:,:)) &
224 /sum(abs(anapot),mask=sim%Mesh%without_ghost_zones%mask3d(:,:,:))
225 sim%aborted=.false.
226 CALL sim%Finalize()
227 print *,"max local relative error: ", pot_err(1)
228 print *,"global relative error: ", pot_err(2)
229 CASE(3)
230 CALL sim%Run()
231 CALL sim%Finalize()
232 END SELECT
233
234 DEALLOCATE(sim)
235 END DO
236
237 IF (ASSOCIATED(anapot)) DEALLOCATE(anapot)
238
239
240CONTAINS
241
242 SUBROUTINE makeconfig(Sim,config)
243 IMPLICIT NONE
244 !------------------------------------------------------------------------!
245 CLASS(fosite), INTENT(INOUT) :: Sim
246 TYPE(dict_typ),POINTER :: config
247 !------------------------------------------------------------------------!
248 ! Local variable declaration
249 TYPE(dict_typ),POINTER :: mesh,boundary,timedisc,datafile,&
250 fluxes,physics,grav,rotframe,sources
251 !------------------------------------------------------------------------!
252 !------------------------------------------------------------------------!
253 physics => dict( &
254 "problem" / euler,&
255 "mu" / mu, &
256 "gamma" / gamma, &
257 "units" / geometrical)
258
259 fluxes => dict( &
260 "fluxtype" / kt, &!HLLC, &
261 "order" / linear, &
262 "variables" / primitive, &
263 "limiter" / monocent, &
264 "theta" / 1.2)
265
266 mesh => dict( &
267 "meshtype" / midpoint, &
268 "geometry" / cylindrical, &
269 "xmin" / rmin, &
270 "xmax" / rmax, &
271! "geometry" / LOGPOLAR, &
272! "xmin" / LOG(RMIN/1.), &
273! "xmax" / LOG(RMAX/1.), &
274 "gparam" / 1., &
275 "omega" / omega, &
276 "inum" / xres, &
277 "jnum" / yres, &
278 "knum" / 1, &
279 "ymin" / (0.), &
280 "ymax" / (2.*pi), &
281 "zmin" / (0.), &
282 "zmax" / (0.), &
283 "decomposition" / (/ -1, 1/))!, &
284! "output/dl" / 1, &
285! "output/radius" / 1, &
286! "output/volume" / 1 )
287
288 boundary => dict( &
289 "western" / noslip, &
290 "eastern" / noslip, &
291 "southern" / periodic, &
292 "northern" / periodic, &
293 "bottomer" / reflecting, &
294 "topper" / reflecting)
295
296 grav => dict( &
297 "stype" / gravity, &
298 "output/accel" / 1, &
299 "self/gtype" / spectral, &
300 "self/green" / 1, &
301 "self/sigma" / sigma)
302! "self/output/potential" / 1)
303
304 rotframe => dict( &
305 "stype" / rotating_frame, &
306 "x" / 0.0, &
307 "y" / 0.0, &
308 "z" / 0.0)
309
310 sources => dict( &
311! "rotframe" / rotframe, &
312 "grav" / grav)
313
314 timedisc => dict( &
315 "method" / ssprk, &
316 "rhstype" / 1, &
317! "dtmax" / 8., &
318! "fargo" / 1, &
319! "method" / RK_FEHLBERG, &
320! "order" / 5, &
321 "tol_rel" / 1.0e-3, &
322 "cfl" / 0.3, &
323 "stoptime" / p0, &
324 "dtlimit" / 1.0e-10, &
325 "maxiter" / 2000000000)
326
327 datafile => dict( &
328 "fileformat" / xdmf, &
329 "filename" / "orbitingcylinders", &
330 "count" / onum)
331
332 config => dict( &
333 "physics" / physics, &
334 "fluxes" / fluxes, &
335 "mesh" / mesh, &
336 "boundary" / boundary, &
337 "sources" / sources, &
338 "timedisc" / timedisc, &
339 "datafile" / datafile)
340
341 END SUBROUTINE makeconfig
342
343
344 SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes,Sources)
345#ifdef HAVE_FGSL
346 USE fgsl
347#endif
348 USE functions, ONLY : ei,bessel_i0,bessel_i1,bessel_k0,bessel_k1,bessel_k0e
349 USE sources_base_mod, ONLY : sources_base
350 USE sources_gravity_mod, ONLY : sources_gravity
351 IMPLICIT NONE
352 !------------------------------------------------------------------------!
353 CLASS(timedisc_base), INTENT(INOUT) :: Timedisc
354 CLASS(mesh_base), INTENT(IN) :: Mesh
355 CLASS(physics_base), INTENT(INOUT) :: Physics
356 CLASS(fluxes_base), INTENT(IN) :: Fluxes
357 CLASS(sources_list), ALLOCATABLE, INTENT(IN) :: Sources
358 !------------------------------------------------------------------------!
359 ! Local variable declaration
360 INTEGER :: i,j,k,dir,ig
361 REAL,DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: r1,r2,r3
362 REAL, DIMENSION(:,:,:), POINTER :: r
363 CLASS(sources_base), POINTER :: sp => null()
364 CLASS(gravity_base), POINTER :: gp => null()
365 !------------------------------------------------------------------------!
366 r => mesh%radius%bcenter
367 ! get gravity spectral solver
368 IF (ALLOCATED(sources)) THEN
369 sp => sources%GetSourcesPointer(gravity)
370 ELSE
371 CALL physics%Error("orbitingcylinders::InitData","no source terms initialized")
372 END IF
373 IF (ASSOCIATED(sp)) THEN
374 SELECT TYPE (sp)
375 CLASS IS(sources_gravity)
376 numpot => sp%pot%data4d(:,:,:,1)
377 gp => sp%glist%GetGravityPointer(spectral)
378 END SELECT
379 ELSE
380 CALL physics%Error("orbitingcylinders::InitData","gravity module not initialized")
381 END IF
382 IF (.NOT.ASSOCIATED(numpot)) &
383 CALL physics%Error("orbitingcylinders::InitData","no pointer to numerical gravitational potential found")
384 IF (.NOT.ASSOCIATED(gp)) &
385 CALL physics%Error("orbitingcylinders::InitData","no spectral gravity solver initialized")
386
387 IF (.NOT.ASSOCIATED(anapot)) ALLOCATE(anapot(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX))
388
389 SELECT CASE(green)
390 CASE(1)
391 SELECT TYPE (pvar => timedisc%pvar)
392 CLASS IS(statevector_euler)
393 r1(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x1(1),x1(2))
394 r2(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x3(1),x3(2))
395 r3(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x4(1),x4(2))
396 pvar%density%data3d(:,:,:) &
397 = 2.0 / (2.*pi*sigma**2) * exp(-sqrt(r1(:,:,:))/sigma) &
398 + 0.5 / (2.*pi*sigma**2) * exp(-sqrt(r2(:,:,:))/sigma) &
399 + 1.0 / (2.*pi*sigma**2) * exp(-sqrt(r3(:,:,:))/sigma)
400 pvar%pressure%data1d(:) = 1.0
401 pvar%velocity%data1d(:) = 0.0
402
403 IF (ASSOCIATED(anapot)) THEN
404 r1(:,:,:) = sqrt(rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x1(1),x1(2))) / (2.*sigma)
405 r2(:,:,:) = sqrt(rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x3(1),x3(2))) / (2.*sigma)
406 r3(:,:,:) = sqrt(rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x4(1),x4(2))) / (2.*sigma)
407 DO k=mesh%KGMIN,mesh%KGMAX
408 DO j=mesh%JGMIN,mesh%JGMAX
409 DO i=mesh%IGMIN,mesh%IGMAX
410 anapot(i,j,k) &
411#ifdef HAVE_FGSL
412 = -2.0 * r1(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r1(i,j,k))*fgsl_sf_bessel_kc1(r1(i,j,k))&
413 -fgsl_sf_bessel_ic1(r1(i,j,k))*fgsl_sf_bessel_kc0(r1(i,j,k))) &
414 -0.5 * r2(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r2(i,j,k))*fgsl_sf_bessel_kc1(r2(i,j,k))&
415 -fgsl_sf_bessel_ic1(r2(i,j,k))*fgsl_sf_bessel_kc0(r2(i,j,k))) &
416 -1.0 * r3(i,j,k)/sigma*(fgsl_sf_bessel_ic0(r3(i,j,k))*fgsl_sf_bessel_kc1(r3(i,j,k))&
417 -fgsl_sf_bessel_ic1(r3(i,j,k))*fgsl_sf_bessel_kc0(r3(i,j,k)))
418#else
419 = -2.0 * r1(i,j,k)/sigma*(bessel_i0(r1(i,j,k))*bessel_k1(r1(i,j,k))&
420 -bessel_i1(r1(i,j,k))*bessel_k0(r1(i,j,k))) &
421 -0.5 * r2(i,j,k)/sigma*(bessel_i0(r2(i,j,k))*bessel_k1(r2(i,j,k))&
422 -bessel_i1(r2(i,j,k))*bessel_k0(r2(i,j,k))) &
423 -1.0 * r3(i,j,k)/sigma*(bessel_i0(r3(i,j,k))*bessel_k1(r3(i,j,k))&
424 -bessel_i1(r3(i,j,k))*bessel_k0(r3(i,j,k)))
425#endif
426 END DO
427 END DO
428 END DO
429 END IF
430 END SELECT
431 SELECT TYPE(sp)
432 CLASS IS(sources_gravity)
433 CALL sp%UpdateGravity(mesh,physics,fluxes,timedisc%pvar,0.0,0.0)
434 END SELECT
435 CASE(2)
436 SELECT TYPE (pvar => timedisc%pvar)
437 CLASS IS(statevector_euler)
438 r1(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x1(1),x1(2))
439 r2(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x3(1),x3(2))
440 r3(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x4(1),x4(2))
441 pvar%density%data3d(:,:,:) &
442 = 2.0 / (2.*pi*sigma**2) * exp(-r1(:,:,:)/(2.*sigma**2)) &
443 + 0.5 / (2.*pi*sigma**2) * exp(-r2(:,:,:)/(2.*sigma**2)) &
444 + 1.0 / (2.*pi*sigma**2) * exp(-r3(:,:,:)/(2.*sigma**2))
445 pvar%pressure%data1d(:) = 1.0
446 pvar%velocity%data1d(:) = 0.0
447
448 IF (ASSOCIATED(anapot)) THEN
449 r1(:,:,:) = rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x1(1),x1(2))
450 r2(:,:,:) = rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x3(1),x3(2))
451 r3(:,:,:) = rsq1(mesh%radius%faces(:,:,:,1),mesh%curv%faces(:,:,:,1,2),x4(1),x4(2))
452 anapot(:,:,:) &
453 = - 2.0 / sqrt(r1(:,:,:)) * erf(sqrt(r1(:,:,:) / 2.) / sigma) &
454 - 0.5 / sqrt(r2(:,:,:)) * erf(sqrt(r2(:,:,:) / 2.) / sigma) &
455 - 1.0 / sqrt(r3(:,:,:)) * erf(sqrt(r3(:,:,:) / 2.) / sigma)
456 END IF
457 END SELECT
458 SELECT TYPE(sp)
459 CLASS IS(sources_gravity)
460 CALL sp%UpdateGravity(mesh,physics,fluxes,timedisc%pvar,0.0,0.0)
461 END SELECT
462 CASE(3)
463 SELECT TYPE (pvar => timedisc%pvar)
464 CLASS IS(statevector_euler)
465 r1(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x1(1),x1(2))
466 r2(:,:,:) = rsq1(r(:,:,:),mesh%curv%bcenter(:,:,:,2),x2(1),x2(2))
467 pvar%density%data3d(:,:,:) = 0.02/(3.2*pi) &
468 + 0.99 * ( exp(-0.5*(r1(:,:,:)/sigma**2))/(2.*pi*sigma**2) &
469 + exp(-0.5*(r2(:,:,:)/sigma**2))/(2.*pi*sigma**2))
470
471 pvar%pressure%data3d(:,:,:) = 0.02/(3.2*pi) &
472 + gn / (2.*pi*sigma**2) &
473 * ( ei(-(r1(:,:,:)/sigma**2)) - ei(-0.5*(r1(:,:,:)/sigma**2)) &
474 + ei(-(r2(:,:,:)/sigma**2)) - ei(-0.5*(r2(:,:,:)/sigma**2)))
475
476 pvar%velocity%data4d(:,:,:,1) = 0.0
477 pvar%velocity%data4d(:,:,:,2) = r(:,:,:)*(1.0 - omega)
478 END SELECT
479 END SELECT
480
481 DO dir=west,east
482 SELECT TYPE(bound => timedisc%Boundary%boundary(dir)%p)
483 CLASS IS (boundary_noslip)
484 DO k=mesh%KMIN,mesh%KMAX
485 DO j=mesh%JMIN,mesh%JMAX
486 DO ig=1,mesh%GNUM
487 SELECT CASE (dir)
488 CASE(west)
489 i = mesh%IMIN-ig
490 CASE(east)
491 i = mesh%IMAX+ig
492 END SELECT
493 bound%data(ig,j,k,physics%YVELOCITY) &
494 = timedisc%pvar%data4d(i,j,k,physics%YVELOCITY)
495 END DO
496 END DO
497 END DO
498 END SELECT
499 END DO
500
501 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
502 END SUBROUTINE initdata
503
504 FUNCTION rsq(Mesh,r,x) RESULT(res)
505 IMPLICIT NONE
506 !------------------------------------------------------------------------!
507 CLASS(mesh_base) :: mesh
508 REAL, DIMENSION(2) :: x
509 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: r
510 REAL, DIMENSION(Mesh%IGMIN:Mesh%IGMAX,Mesh%JGMIN:Mesh%JGMAX,Mesh%KGMIN:Mesh%KGMAX) :: res
511 !------------------------------------------------------------------------!
512 INTENT(IN) :: x
513 !------------------------------------------------------------------------!
514 res = r**2 + x(1)**2 - 2.*r*x(1) * cos(mesh%curv%bcenter(:,:,:,2)-x(2))
515 END FUNCTION rsq
516
517 ELEMENTAL FUNCTION rsq1(r,phi,r0,phi0) RESULT(res)
518 IMPLICIT NONE
519 !------------------------------------------------------------------------!
520 REAL, INTENT(IN) :: r, phi, r0, phi0
521 REAL :: res
522 !------------------------------------------------------------------------!
523 res = r**2 + r0**2 - 2.*r*r0 * cos(phi-phi0)
524 END FUNCTION rsq1
525
526END PROGRAM orbitingcylinders
initdata
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
makeconfig
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
fosite_mod
Definition: fosite.f90:41
functions
Mathematical functions.
Definition: functions.f90:48
functions::bessel_k0e
elemental real function, public bessel_k0e(x)
Compute the exponential scaled modified Bessel function of the second kind e.g. K0e = EXP(x) * K0(x) ...
Definition: functions.f90:692
functions::bessel_i1
elemental real function, public bessel_i1(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
Definition: functions.f90:615
functions::ei
elemental real function, public ei(x)
Computes the exponential integral Ei(x) for x != 0 Ei(x) := int_-oo^x exp(t)/t dt Implementation simi...
Definition: functions.f90:758
functions::bessel_k1
elemental real function, public bessel_k1(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
Definition: functions.f90:721
functions::bessel_k0
elemental real function, public bessel_k0(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
Definition: functions.f90:657
functions::pi
real, parameter pi
Definition: functions.f90:52
functions::bessel_i0
elemental real function, public bessel_i0(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
Definition: functions.f90:577
sources_base_mod
generic source terms module providing functionaly common to all source terms
Definition: sources_base.f90:42
sources_gravity_mod
generic gravity terms module providing functionaly common to all gravity terms
Definition: sources_gravity.f90:41
rsq1
elemental real function rsq1(r, phi, r0, phi0)
Definition: orbitingcylinders.f90:518
orbitingcylinders
program orbitingcylinders
Definition: orbitingcylinders.f90:153
rsq
real function, dimension(mesh%igmin:mesh%igmax, mesh%jgmin:mesh%jgmax, mesh%kgmin:mesh%kgmax) rsq(Mesh, r, x)
Definition: orbitingcylinders.f90:505
omega
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
Definition: sources_viscosity.f90:432
fosite_mod::fosite
main fosite class
Definition: fosite.f90:71
sources_base_mod::sources_base
Definition: sources_base.f90:57
sources_gravity_mod::sources_gravity
Definition: sources_gravity.f90:58
Generated by   doxygen 1.9.4