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. /)
165 REAL,
PARAMETER :: gn = 1.
166 REAL,
PARAMETER :: p0 = 1.
167 REAL,
PARAMETER :: p = 10.
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.
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
181 CLASS(
fosite),
ALLOCATABLE :: sim
183 CHARACTER(LEN=1) :: str
184 REAL,
DIMENSION(2) :: pot_err
185 REAL,
DIMENSION(:,:,:),
POINTER :: numpot => null(), anapot => null()
190 CALL sim%InitFosite()
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)
201 CALL setattr(sim%config,
"/timedisc/stoptime",1e-5*p0)
202 CALL setattr(sim%config,
"/datafile/count",1)
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)
213 CALL initdata(sim%Timedisc,sim%Mesh,sim%Physics,sim%Fluxes,sim%Sources)
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(:,:,:))
227 print *,
"max local relative error: ", pot_err(1)
228 print *,
"global relative error: ", pot_err(2)
237 IF (
ASSOCIATED(anapot))
DEALLOCATE(anapot)
245 CLASS(
fosite),
INTENT(INOUT) :: Sim
246 TYPE(dict_typ),
POINTER :: config
249 TYPE(dict_typ),
POINTER :: mesh,boundary,timedisc,datafile,&
250 fluxes,physics,grav,rotframe,sources
257 "units" / geometrical)
262 "variables" / primitive, &
263 "limiter" / monocent, &
267 "meshtype" / midpoint, &
268 "geometry" / cylindrical, &
283 "decomposition" / (/ -1, 1/))
289 "western" / noslip, &
290 "eastern" / noslip, &
291 "southern" / periodic, &
292 "northern" / periodic, &
293 "bottomer" / reflecting, &
294 "topper" / reflecting)
298 "output/accel" / 1, &
299 "self/gtype" / spectral, &
301 "self/sigma" / sigma)
305 "stype" / rotating_frame, &
321 "tol_rel" / 1.0e-3, &
324 "dtlimit" / 1.0e-10, &
325 "maxiter" / 2000000000)
328 "fileformat" / xdmf, &
329 "filename" /
"orbitingcylinders", &
333 "physics" / physics, &
336 "boundary" / boundary, &
337 "sources" / sources, &
338 "timedisc" / timedisc, &
339 "datafile" / datafile)
344 SUBROUTINE initdata(Timedisc,Mesh,Physics,Fluxes,Sources)
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
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
364 CLASS(gravity_base),
POINTER :: gp => null()
366 r => mesh%radius%bcenter
368 IF (
ALLOCATED(sources))
THEN
369 sp => sources%GetSourcesPointer(gravity)
371 CALL physics%Error(
"orbitingcylinders::InitData",
"no source terms initialized")
373 IF (
ASSOCIATED(sp))
THEN
376 numpot => sp%pot%data4d(:,:,:,1)
377 gp => sp%glist%GetGravityPointer(spectral)
380 CALL physics%Error(
"orbitingcylinders::InitData",
"gravity module not initialized")
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")
387 IF (.NOT.
ASSOCIATED(anapot))
ALLOCATE(anapot(mesh%IGMIN:mesh%IGMAX,mesh%JGMIN:mesh%JGMAX,mesh%KGMIN:mesh%KGMAX))
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
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
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)))
433 CALL sp%UpdateGravity(mesh,physics,fluxes,timedisc%pvar,0.0,0.0)
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
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))
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)
460 CALL sp%UpdateGravity(mesh,physics,fluxes,timedisc%pvar,0.0,0.0)
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))
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)))
476 pvar%velocity%data4d(:,:,:,1) = 0.0
477 pvar%velocity%data4d(:,:,:,2) = r(:,:,:)*(1.0 -
omega)
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
493 bound%data(ig,j,k,physics%YVELOCITY) &
494 = timedisc%pvar%data4d(i,j,k,physics%YVELOCITY)
501 CALL physics%Convert2Conservative(timedisc%pvar,timedisc%cvar)
504 FUNCTION rsq(Mesh,r,x)
RESULT(res)
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
514 res = r**2 + x(1)**2 - 2.*r*x(1) * cos(mesh%curv%bcenter(:,:,:,2)-x(2))
517 ELEMENTAL FUNCTION rsq1(r,phi,r0,phi0)
RESULT(res)
520 REAL,
INTENT(IN) :: r, phi, r0, phi0
523 res = r**2 + r0**2 - 2.*r*r0 * cos(phi-phi0)
subroutine initdata(Mesh, Physics, Fluxes, Timedisc)
subroutine makeconfig(Sim, config)
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) ...
elemental real function, public bessel_i1(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
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...
elemental real function, public bessel_k1(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
elemental real function, public bessel_k0(x)
Compute the modified Bessel function of the second kind as polynomial expansion, which has been propo...
elemental real function, public bessel_i0(x)
Compute the modified Bessel function of the first kind as polynomial expansion, which has been propos...
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)