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!----------------------------------------------------------------------------!
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
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
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
Definition: bondi2d.f90:274
subroutine makeconfig(Sim, config)
Definition: bondi2d.f90:165
Mathematical functions.
Definition: functions.f90:48
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
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
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
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
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
real, parameter pi
Definition: functions.f90:52
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
generic source terms module providing functionaly common to all source terms
generic gravity terms module providing functionaly common to all gravity terms
elemental real function rsq1(r, phi, r0, phi0)
program orbitingcylinders
real function, dimension(mesh%igmin:mesh%igmax, mesh%jgmin:mesh%jgmax, mesh%kgmin:mesh%kgmax) rsq(Mesh, r, x)
pure real function, dimension(size(ephir%data2d, dim=1)) omega(ephir, velocity)
main fosite class
Definition: fosite.f90:71