adam_flume_ic_object

ADAM, FLUME initial conditions class definition.

Accepted [initial_conditions].(type):

• uniform: the state of region 1 everywhere, density and pressure perturbed by the signed relative amplitude s, rho (1 + s h1), p (1 + s h2), with h1, h2 in [-1, 1) hashed from the block Morton code and the cell indexes: the same on every rank decomposition and on every backend (no random generator state);

• isentropic-vortex: the free stream of region 1 plus an isentropic vortex of axis z centred at (x0, y0), of radius radius and velocity strength strength (Shu 1998, ICASE 97-65, section 5.1, written for any free stream); optional shock_x: past it the Rankine-Hugoniot state of a stationary shock, region 1 upstream and supersonic: du = -strength/(2 pi) (y-y0)/radius exp((1-r^2)/2), dv = strength/(2 pi) (x-x0)/radius exp((1-r^2)/2), d(p/rho) = -(gamma-1)/gamma strength^2/(8 pi^2) exp(1-r^2), isentropic, r = |x-x0| / radius; an exact steady solution convected by the free stream;

• riemann-problem: piecewise-constant axis-aligned regions, a cell belongs to a region iff emin < center <= emax on every axis. A cell covered by no region is fatal: leaving it at zero density would divide by zero downstream;

• glm-pulse (MHD only, the GLM d'Alembert test, issue #41, MV-3): the state of region 1 plus a Gaussian pulse in the magnetic field along pulse_axis (x, y or z), B_axis += pulse_amplitude exp(-((s - pulse_center) / pulse_width)^2), s the cell coordinate along the axis; pressure and velocity unchanged (psi zero);

• divb-peak (MHD only, the Dedner et al. 2002 peak in B_x, issue #41, MV-10): the state of region 1 plus B_x += peak_amplitude (1 - (s / peak_radius)^2)^2 for s < peak_radius, s the distance in the x-y plane from (peak_x0, peak_y0) (Dedner: peak_amplitude = 1 / sqrt(4 pi), peak_radius = 1/8); a non-zero initial div(B);

• mhd-linear-wave (MHD only, issue #41, MV-5; Stone et al. 2008, section 8.2): the state q0 of region 1 (global frame) plus wave_amplitude R sin(2 pi (x cos a + y sin a) / wavelength), a = wave_angle (degrees, in the x-y plane), R the right eigenvector of the right-going wave (fast, alfven, slow or entropy) of the conservative system at q0 in the frame of the wave normal (the library Roe-Balsara core, mhd_eigenvectors), rotated back to x, y, z; after one period wavelength / |u_n + c| the exact solution is the initial state;

• mhd-cpaw (MHD only, issue #41, MV-6; Toth 2000, J. Comput. Phys. 161): the circularly polarised Alfven wave, an exact nonlinear solution of any amplitude. Density and pressure of region 1 (whose velocity and field must be zero), the field b_par n + wave_amplitude (sin(phi) t1 + h cos(phi) t2) and the velocity -sign(b_par) B_perp / sqrt(rho), phi = 2 pi (x cos a + y sin a) / wavelength, a = wave_angle (degrees), n = (cos a, sin a, 0), t1 = (-sin a, cos a, 0), t2 = z, h = +1 (polarisation = right) or -1 (left, the mirror image under z -> -z); |B_perp| is uniform, so the total pressure is too, and the wave travels along +n unchanged at |b_par| / sqrt(rho): after one period wavelength sqrt(rho) / |b_par| the exact solution is the initial state;

• mhd-vortex (MHD only, issue #41, MV-7; Balsara 2004, ApJS 151): the free stream of region 1 (bx = by = 0, fatal otherwise) plus a magnetised vortex of axis z centred at (x0, y0), of radius radius: dv = kappa/(2 pi) e (-(y-y0), x-x0) / radius, dB = mu/(2 pi) e (-(y-y0), x-x0) / radius, dp = (mu^2 (1-r^2) - rho kappa^2) / (8 pi^2) e^2, e = exp((1-r^2)/2), r = |x-x0| / radius, density uniform: the radial balance dp/dr = rho v^2/r - B^2/r - d(B^2/2)/dr (rationalised units) holds exactly, so the vortex is a steady solution convected by the free stream, with a divergence-free field;

• orszag-tang (MHD only, issue #41, MV-12; Stone et al. 2008, section 8.4): the Orszag-Tang vortex on the unit period, keyless and with no region sections: rho = 25/(36 pi), p = 5/(12 pi), u = -sin(2 pi y), v = sin(2 pi x), bx = -sin(2 pi y) / sqrt(4 pi), by = sin(4 pi x) / sqrt(4 pi), w = bz = 0 (divergence-free pointwise, symmetric under the 180 degrees rotation about (1/2, 1/2)). Evaluated as odd functions of s = 2 x - 1, t = 2 y - 1 (u = sin(pi t), v = -sin(pi s), by = sin(2 pi s) / sqrt(4 pi)): on cell centres that are exact binary fractions the rotated cell has exactly -s, -t, so the initial state is bitwise symmetric;

• mhd-rotor (MHD only, issue #41, MV-13; Balsara and Spicer 1999, Toth 2000 first rotor): the ambient state of region 1 plus a dense disc centred at (x0, y0) in rigid rotation: for r < r0 density rho_in and velocity v0 (-(y-y0), x-x0) / r0 (speed v0 at r0); for r0 <= r < r1 the linear taper f = (r1 - r) / (r1 - r0), density rho + (rho_in - rho) f, velocity f v0 (-(y-y0), x-x0) / r; pressure and field of region 1 everywhere; the velocities add to the ambient one;

• field-loop (MHD only, issue #41, MV-9; Gardiner and Stone 2005, Mignone and Tzeferacos 2010 section 4.4.1): the state of region 1 plus the field of the vector potential A_z = loop_amplitude (loop_radius - r) for r < loop_radius, r = |x - (x0, y0)| in the x-y plane: dB = loop_amplitude (-(y-y0), x-x0) / r inside the loop, zero outside (divergence-free pointwise, discontinuous at r = loop_radius); with a uniform w the discrete div(B) feeds B_z (d B_z / dt = w div(B));

• rotated-riemann (issue #41, MV-8; Toth 2000, J. Comput. Phys. 161, section 6.3.2): a one-dimensional Riemann problem rotated in the x-y plane, periodic along its normal. The normal is (normal_x, normal_y), the coordinate s = normal_x x + normal_y y (unnormalised), and with f = modulo(s - interface_1, period), d = interface_2 - interface_1 and w = interface_2_width a cell is in region 2 iff 0 < f <= d, in the linear ramp of the primitive variables from region 2 to region 1 iff d < f <= d + w (at f = d + w region 1), in region 1 otherwise: a jump at interface_1, and at interface_2 a jump (w = 0) or a ramp. On a periodic line the velocity jump of a shock tube must be undone somewhere: a second jump is the mirror problem, an expansion that may approach vacuum (for Ryu-Jones 1a: density ~2e-4, Alfven speed ~100); the ramp spreads it. The region states are given in the frame of the normal (u, v = normal and tangential velocity, bx, by = normal and tangential field, w, bz along z) and rotated to x, y. With integer normal components on a square periodic domain of side period, the state is periodic in x and in y (Toth's shifted-periodic strip is not needed).

• channel (issue #119, T3.1; euler): a plane channel between y = +-ch_h, ch_rho, bulk velocity ch_ub, wall Mach ch_mach, ch_lx, ch_lz, perturbation ch_pert.

• boundary-layer (issue #119, T3.7; euler): the mean profile of a table (the format of the digital-filter inflow, adam_flume_inflow_df_object) over a wall at y = bl_wall (default 0), flow along +x: bl_profile, bl_delta, free stream bl_rho, bl_u, bl_p; uniform along x and z.

• taylor-green (issue #119, T0.2; euler, euler-de): the compressible Taylor-Green vortex of DeBonis (2013), tg_density (rho_0), tg_velocity (V_0), tg_mach (M_0), tg_length (L, default 1), no regions: u = V_0 sin(x/L) cos(y/L) cos(z/L), v = -V_0 cos(x/L) sin(y/L) cos(z/L), w = 0, p = p_0 + rho_0 V_0^2 / 16 (cos(2x/L) + cos(2y/L)) (cos(2z/L) + 2), p_0 = rho_0 V_0^2 / (gamma M_0^2), at uniform temperature, rho = rho_0 p / p_0. The box is [0, 2 pi L]^3, periodic; the initial condition does not check it.

• isotropic-turbulence (issue #119, T0.5; euler, euler-de): the decaying compressible isotropic turbulence of Johnsen et al. (2010), hit_density (rho_0), hit_urms (u_rms), hit_mach (M_t = sqrt(3) u_rms / c_0), hit_k0 (k_0), hit_seed (default 1), no regions: uniform density and pressure, p_0 = rho_0 c_0^2 / gamma, and a solenoidal velocity u = sum 2 Re(a_k exp(i k.x)) over the integer wavevectors of a half space with |k| <= 3 k_0, of random phases (Rogallo 1981), each shell s = nint(|k|) holding E(s) = 16 sqrt(2/pi) u_rms^2 / k_0 (s/k_0)^4 exp(-2 (s/k_0)^2), the sum rescaled so that the mean of |u|^2 is exactly 3 u_rms^2. The sequence depends on the seed alone (MINSTD, its own state), so the field is the same on any rank count, block layout and resolution. The box is [0, 2 pi]^3, periodic; the initial condition does not check it, nor that 3 k_0 is below the Nyquist wavenumber of the grid.

• shu-osher (Euler only, issue #47, RV-3; Shu & Osher 1989, J. Comput. Phys. 83, example 8): a shock running into a density wave along axis (x, y or z), coordinate s: the state of region 1 for s <= interface, else the state of region 2 with density r + rho_amplitude sin(rho_wavenumber s) (the classic problem: region 1 3.857143, 2.629369, 10.33333, region 2 1, 0, 1, interface -4, amplitude 0.2, wavenumber 5 on [-5, 5], t = 1.8). The regions take no extent keys;

• linear (verification, issue #65 P0): every conservative variable linear in space, q = q_1 (1 + g . x) with q_1 the state of region 1 and the relative gradient g = (gradient_x, gradient_y, gradient_z). The same-level copy, the 2:1 restriction (mean of the fine cells) and the tricubic coarse->fine fill reproduce it exactly, so the ghost cells of a step-0 checkpoint can be checked against it (verification/ghosts/ghost_probe.py). The state is admissible only while 1 + g . x > 0 over the domain; the probe runs no step.

• polynomial (verification, issue #68 D1): q = q_1 (1 + sum_d g_d x_d^n), n = degree (an integer >= 1; n = 1 is linear). Per-direction degree n: exact under the same-level copy and the boundary conditions (copies of cells), under the 2:1 restriction when it is exact on degree n (lagrange6: n <= 5; the mean: n = 1) and under the coarse->fine fill likewise (tricubic: n <= 3, quintic: n <= 5), so ghost_probe.py checks the high-order seam operators in situ. Not available under a non-dimensional reference.

The primitive keys of a region follow the physical model: r, u, v, w, p (Euler), plus bx, by, bz (MHD; psi is zero). isentropic-vortex is Euler only (issue #41, section 3.7).

• isothermal-atmosphere, polytropic-atmosphere (issue #73, M7-G0): the hydrostatic equilibrium of the constant gravity g of [gravity] (required), with the state of region 1 at the origin (its velocity and field uniform): isothermal, rho = rho_1 exp(g . x / T_1), p = rho T_1, T_1 = p_1 / rho_1; polytropic, p = p_1 (rho / rho_1)^G with G = polytropic_gamma > 1, rho = rho_1 (1 + (G - 1) / G rho_1 / p_1 g . x)^(1 / (G - 1)), fatal where the bracket is not positive (the atmosphere's top lies inside the domain). Point values at the cell centres.
• rayleigh-taylor (issue #73, M7-G3, GV-4): two isothermal layers under the constant gravity g of [gravity] (required, along one axis), cold and dense over hot and light, joined by a tanh temperature at the height rt_height (h = -g . x / |g|): T = a + b tanh((h - h_0) / rt_width), a = (T_l + T_h) / 2, b = (T_h - T_l) / 2, T_l = p_1 / rho_1 (region 1, the light layer), T_h = T_l / rt_density_ratio; hydrostatic, p = p_1 exp(-|g| int_{h_0}^h dh / T) in closed form (p(h_0) = p_1), rho = p / T, at rest with region 1's field. Seeded with a solenoidal single mode, v_h = A cos(k x_t) f(h) along -g and u_t = -A / k sin(k x_t) f'(h) along x_t (the first axis other than g's), A = rt_amplitude, k = 2 pi / rt_wavelength, f = exp(-((h - h_0) / rt_seed_width)^2): divergence-free, so the seed starts no sound waves. Every layer is stably stratified, so the interface is the only instability; the linear growth rate is verification/gravity/rt_oracle.py.
• polytrope (issue #81, M8-P4, SV-7/SV-8; Euler, [gravity].(self_gravity) required): an n = 1 polytrope star of radius star_radius (a) and central density star_density (rho_c), centred at (star_x0, star_y0, star_z0) (default the origin), in a light atmosphere set by region 1: inside, rho = rho_1 + rho_s, p = p_1 + f K rho_s^2, rho_s = rho_c sin(k r) / (k r) (k = pi / a), K = 2 G a^2 / pi (the star alone is then in equilibrium: Lane-Emden n = 1) and f = star_pressure_factor (default 1; f < 1 collapses it); outside, the atmosphere is isothermal at T = p_1 / rho_1 and hydrostatic in the star's field, rho = rho_1 exp(G M / T (1/r - 1/a)), p = rho T, M = 4 rho_c a^3 / pi (a uniform atmosphere would fall onto the star and thin at its surface, which only the positivity limiter survives); region 1's velocity. Density and pressure are continuous at the surface.
• singular-isothermal-sphere (issue #94, M12-P4, SI si11; Euler, [gravity].(self_gravity) required): the singular isothermal sphere at rest centred at (sis_x0, sis_y0, sis_z0) (default the origin), rho = A a^2 / (4 pi G r^2) (the cell mean), p = rho a^2, a^2 = p_1 / rho_1 (region 1's temperature: only the ratio is used), region 1's velocity, A = sis_overdensity (default 2, > 0): A = 2 is the hydrostatic sphere of Shu (1977), whose inside-out collapse accretes at 0.975 a^3 / G; A > 2 collapses everywhere at m_0(A) a^3 / G (Shu's A-family; A = 29.3 is the run of Federrath et al. 2010, section 3.4). Fatal at a cell centre on the centre (place it on a cell corner). With [eos] type = isothermal (issue #100, P4) a^2 = c0^2, the gas's own sound speed.
• bonnor-ebert (issue #100, P4; [eos] type = isothermal and [gravity].(self_gravity) required): the isothermal Lane-Emden sphere at rest centred at (be_x0, be_y0, be_z0) (default the origin), rho = rho_c exp(-psi(xi)), xi = r / r_0, r_0 = c0 / sqrt(4 pi G rho_c), rho_c = be_density (> 0), psi'' + 2 psi' / xi = exp(-psi), psi(0) = psi'(0) = 0 (RK4 on xi steps of 1e-3, read by cubic Hermite interpolation), region 1's velocity. It fills the domain: closed walls confine it, and the sphere inscribed in the box is a Bonnor-Ebert sphere of its xi, stable below 6.45 (Bonnor 1956, Ebert 1955). be_overdensity (default 1, > 0) multiplies the density everywhere (above 1 the sphere is no equilibrium and collapses: a supercritical sphere left exactly in equilibrium may as well expand).
• boss-bodenheimer (issue #100, P4; [eos] isothermal or barotropic, [gravity].(self_gravity) required): the rotating cloud of Boss & Bodenheimer (1979) as set by Burkert & Bodenheimer (1993): a sphere of radius bb_radius centred at (bb_x0, bb_y0, bb_z0) (default the origin), rho = bb_density (1 + bb_amplitude cos(2 phi)) (bb_amplitude default 0.1, phi the azimuth about z), in solid-body rotation about z at bb_omega; outside, region 1's density and velocity (the light ambient medium). The density and velocity only: the pressure is the gas's own P(rho).

Source: src/app/flume/common/adam_flume_ic_object.F90

Dependencies

Contents

• flume_ic_object
• initialize
• load_from_file
• set_equilibrium
• set_initial_conditions
• hash_cell
• linear_wave_eigenvector
• linear_wave_eigenvector_inert
• mhd_vortex
• mhd_rotor
• atmosphere_state
• rotate_to_global
• set_isotropic_turbulence
• isotropic_turbulence_modes
• channel_state
• taylor_green
• isentropic_vortex
• becker_shock
• is_barotropic_safe
• description
• lane_emden_isothermal_table
• lane_emden_isothermal_psi
• inverse_square_mean

Variables

Name	Type	Attributes	Description
INI_SECTION_NAME	character(len=18)	parameter	INI section name.
IC_UNIFORM_STR	character(len=7)	parameter	Uniform state.
IC_ISENTROPIC_VORTEX_STR	character(len=17)	parameter	Isentropic vortex in a free stream.
IC_RIEMANN_PROBLEM_STR	character(len=15)	parameter	Piecewise-constant regions.
IC_GLM_PULSE_STR	character(len=9)	parameter	Gaussian pulse in B along an axis (MHD).
PULSE_KEY	character(len=15)	parameter	GLM pulse keys.
IC_DIVB_PEAK_STR	character(len=9)	parameter	Dedner peak in B_x (MHD).
PEAK_KEY	character(len=14)	parameter	div(B) peak keys.
IC_MHD_LINEAR_WAVE_STR	character(len=15)	parameter	MHD linear wave (MHD).
WAVE_KEY	character(len=14)	parameter	Linear wave and CPAW keys.
IC_MHD_CPAW_STR	character(len=8)	parameter	Circularly polarised Alfven wave (MHD).
IC_MHD_VORTEX_STR	character(len=10)	parameter	Magnetised vortex in a free stream (MHD).
IC_ORSZAG_TANG_STR	character(len=11)	parameter	Orszag-Tang vortex (MHD).
IC_MHD_ROTOR_STR	character(len=9)	parameter	MHD rotor (MHD).
IC_FIELD_LOOP_STR	character(len=10)	parameter	Advected field loop (MHD).
LOOP_KEY	character(len=14)	parameter	Field loop keys.
IC_ROTATED_RIEMANN_STR	character(len=15)	parameter	Rotated periodic Riemann problem.
IC_CHANNEL_STR	character(len=7)	parameter	Plane channel, parabolic with a perturbation (issue #119, T3.1).
IC_BOUNDARY_LAYER_STR	character(len=14)	parameter	Mean profile of a boundary layer (issue #119, T3.7).
BL_KEY	character(len=8)	parameter	Boundary-layer keys.
IC_TAYLOR_GREEN_STR	character(len=12)	parameter	Compressible Taylor-Green vortex (issue #119).
CHANNEL_KEY	character(len=7)	parameter	Channel keys.
TG_KEY	character(len=11)	parameter	Taylor-Green required keys.
IC_ISOTROPIC_TURB_STR	character(len=20)	parameter	Decaying isotropic turbulence (issue #119).
HIT_KEY	character(len=11)	parameter	Its required keys.
IC_SHU_OSHER_STR	character(len=9)	parameter	Shock-density wave interaction (Euler).
IC_LINEAR_STR	character(len=6)	parameter	Linear field (verification of ghost fills).
IC_POLYNOMIAL_STR	character(len=10)	parameter	Polynomial field (verification, #68 D1).
IC_SINE_WAVE_STR	character(len=9)	parameter	Small sine wave on a uniform state (issue #65).
SINE_MODE_SHEAR	character(len=5)	parameter	Sine wave: transverse velocity (shear) wave.
SINE_MODE_ACOUSTIC	character(len=8)	parameter	Sine wave: right-running acoustic wave.
SINE_MODE_MAGNETIC	character(len=8)	parameter	Sine wave: transverse field (MHD).
SINE_MODE_ENTROPY	character(len=7)	parameter	Sine wave: density (entropy) wave.
SINE_MODE_JEANS	character(len=5)	parameter	Sine wave: Jeans mode (self-gravity, #81).
ATM_PULSE_KEY	character(len=15)	parameter	Atmosphere pulse keys.
IC_COUETTE_STR	character(len=7)	parameter	Exact compressible Couette flow (issue #65).
COUETTE_KEY	character(len=16)	parameter	Couette keys.
IC_BECKER_STR	character(len=12)	parameter	Becker's exact viscous shock (issue #65).
IC_ISOTHERMAL_ATM_STR	character(len=21)	parameter	Isothermal atmosphere (issue #73).
IC_POLYTROPIC_ATM_STR	character(len=21)	parameter	Polytropic atmosphere (issue #73).
IC_RAYLEIGH_TAYLOR_STR	character(len=15)	parameter	Two isothermal layers, heavy over light (#73).
IC_POLYTROPE_STR	character(len=9)	parameter	n = 1 polytrope star, self-gravity (#81).
IC_SIS_STR	character(len=26)	parameter	Singular isothermal sphere, self-gravity (#94).
SIS_KEY	character(len=15)	parameter	Singular isothermal sphere keys.
IC_BE_STR	character(len=12)	parameter	Isothermal Lane-Emden sphere, self-gravity (#100).
BE_KEY	character(len=14)	parameter	Bonnor-Ebert keys.
BE_DXI	real(kind=R8P)	parameter	Bonnor-Ebert: Lane-Emden table step in xi.
IC_BB_STR	character(len=16)	parameter	Rotating cloud, m = 2 perturbation (#100).
BB_KEY	character(len=12)	parameter	Boss-Bodenheimer keys.
STAR_KEY	character(len=20)	parameter	Polytrope keys.
RT_KEY	character(len=16)	parameter	RT keys.
BECKER_KEY	character(len=8)	parameter	Becker shock keys.
LINEAR_KEY	character(len=10)	parameter	Linear field gradient keys.
SHU_OSHER_KEY	character(len=14)	parameter	Shu-Osher keys.
ROTATED_KEY	character(len=17)	parameter	Rotated Riemann keys.
ROTOR_KEY	character(len=6)	parameter	MHD rotor keys.
MHD_VORTEX_KEY	character(len=6)	parameter	Magnetised vortex keys.
VORTEX_KEY	character(len=8)	parameter	Vortex keys.
PI	real(kind=R8P)	parameter	Pi greek.
PRIM_KEY	character(len=2)	parameter	Region primitive state keys (MHD: all 8).
EXTENT_KEY	character(len=6)	parameter	Region extents keys.

Derived Types

flume_ic_object

FLUME initial conditions class definition.

Components

Name	Type	Attributes	Description
amr_iterations	integer(kind=I4P)		AMR iterations performed while imposing the initial conditions.
ic_type	character(len=:)	allocatable	Initial conditions type.
regions_number	integer(kind=I4P)		Regions number.
model	integer(kind=I4P)		Physical model id.
nprim	integer(kind=I4P)		Primitive keys number of a region: 5 (Euler) or 8 (MHD).
q_region	real(kind=R8P)	allocatable	Conservative state of each region [nv, regions_number].
scalars	real(kind=R8P)	allocatable	Passive scalars of each region [NS_MAX, regions_number] (#92).
passive_scalars	integer(kind=I4P)		Passive scalars in use (issue #92).
iq_scalar1	integer(kind=I4P)		Slot of the first passive scalar.
emin	real(kind=R8P)	allocatable	Minimum corner of each region [3, regions_number].
emax	real(kind=R8P)	allocatable	Maximum corner of each region [3, regions_number].
gamma	real(kind=R8P)		Specific heats ratio.
eos_arg	real(kind=R8P)		EOS parameter of the barotropic models (issue #100): c0^2 or
prim_1	real(kind=R8P)		Primitive state of region 1, the free stream (first nprim used).
s	real(kind=R8P)		Uniform: relative amplitude of the seeded perturbation.
channel	real(kind=R8P)		Channel: density, bulk velocity, wall Mach number, half height, Lx, Lz, perturbation.
bl	real(kind=R8P)		Boundary layer: delta, free-stream density, speed, pressure, wall position.
bl_prof	real(kind=R8P)	allocatable	Boundary layer: profile table [7, rows].
tg	real(kind=R8P)		Taylor-Green: density, velocity, Mach number, length.
hit	real(kind=R8P)		Isotropic turbulence: density, u_rms, turbulent Mach number, k_0.
hit_seed	integer(kind=I4P)		Isotropic turbulence: seed of the random phases.
vortex	real(kind=R8P)		Isentropic vortex: x0, y0, radius, strength.
has_shock	logical		Isentropic vortex: a stationary shock at shock_x (issue #119, T2.7).
shock_x	real(kind=R8P)		Isentropic vortex: position of the stationary shock.
pulse_axis	integer(kind=I4P)		GLM pulse: axis, 1=x, 2=y, 3=z.
pulse	real(kind=R8P)		GLM pulse: center, width, amplitude.
peak	real(kind=R8P)		div(B) peak: x0, y0, radius, amplitude.
wave	character(len=:)	allocatable	Linear wave family: fast, alfven, slow, entropy.
wave_par	real(kind=R8P)		Linear wave, CPAW: angle (degrees), amplitude, wavelength.
wave_env	real(kind=R8P)		Sine wave: Gaussian envelope width (0: none) and centre along n (issue #121).
wave_r	real(kind=R8P)		Linear wave: right eigenvector, global frame.
b_par	real(kind=R8P)		CPAW: field along the wave normal.
polarisation	real(kind=R8P)		CPAW: +1 right, -1 left.
mvortex	real(kind=R8P)		Magnetised vortex: x0, y0, radius, kappa, mu.
rotor	real(kind=R8P)		MHD rotor: x0, y0, r0, r1, rho_in, v0.
loop	real(kind=R8P)		Field loop: x0, y0, radius, amplitude.
rotated	real(kind=R8P)		Rotated Riemann: normal x, y, interface 1, 2, period, width 2.
prim_2	real(kind=R8P)		Primitive state of region 2 (global frame, first nprim used).
shu_osher_axis	integer(kind=I4P)		Shu-Osher: axis, 1=x, 2=y, 3=z.
shu_osher	real(kind=R8P)		Shu-Osher: interface, density amplitude, density wavenumber.
gradient	real(kind=R8P)		Linear: relative gradient g, q = q_1 (1 + g . x).
degree	integer(kind=I4P)		Polynomial: degree n, q = q_1 (1 + sum_d g_d x_d^n).
sine_mode	character(len=:)	allocatable	Sine wave: shear, acoustic, entropy, magnetic or jeans.
four_pi_g	real(kind=R8P)		Sine wave, jeans: 4 pi G of the self-gravity.
couette	real(kind=R8P)		Couette: fixed wall y, height, moving wall velocity, T_wall.
couette_heat	real(kind=R8P)		Couette: mu U^2 / (k H), the temperature rise scale.
R	real(kind=R8P)		Gas constant.
becker	real(kind=R8P)		Becker shock: upstream Mach number, shock position.
becker_mu	real(kind=R8P)		Becker shock: viscosity.
cp	real(kind=R8P)		Specific heat at constant pressure.
g	real(kind=R8P)		Atmospheres: gravity.
polytropic_gamma	real(kind=R8P)		Polytropic atmosphere: exponent G of p ~ rho^G.
atm_pulse	real(kind=R8P)		Atmospheres: pressure pulse amplitude (relative), centre, width.
rt	real(kind=R8P)		seed amplitude, seed width.
rt_axis	integer(kind=I4P)		Rayleigh-Taylor: axis of the gravity.
rt_t1	integer(kind=I4P)		Rayleigh-Taylor: axis of the seeded mode.
star	real(kind=R8P)		Polytrope: radius, central density, centre, pressure factor.
star_k	real(kind=R8P)		Polytrope: K of p = K rho^2 in equilibrium, 2 G a^2 / pi.
star_gm	real(kind=R8P)		Polytrope: G M of the star, 4 G rho_c a^3 / pi.
sis	real(kind=R8P)		Singular isothermal sphere: centre, overdensity A, G.
be	real(kind=R8P)		Bonnor-Ebert sphere: centre, central density, overdensity,
bb	real(kind=R8P)		Boss-Bodenheimer cloud: radius, density, omega, amplitude,

Type-Bound Procedures

Name	Attributes	Description
description	pass(self)	Return pretty-printed object description.
initialize	pass(self)	Initialize initial conditions.
load_from_file	pass(self)	Load config from file.
set_initial_conditions	pass(self)	Set initial conditions on the blocks interior.
set_equilibrium	pass(self)	Set the hydrostatic equilibrium on the blocks interior.
is_barotropic_safe	pass(self)	True if the IC type builds its density without gamma (#100).
atmosphere_state	pass(self)	Primitive state of the atmosphere at a point.
set_isotropic_turbulence	pass(self)	Set the isotropic turbulence field (#119).

Subroutines

initialize

Initialize initial conditions.

subroutine initialize(self, file_parameters, physics)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	inout		Initial conditions.
file_parameters	type(file_ini)	in		Simulation parameters ini file handler.
physics	type(flume_physics_object)	in		Physics (for the regions state conversion).

Call graph

load_from_file

Load config from file; every key used by the selected type is required and an unknown type is fatal.

subroutine load_from_file(self, file_parameters, physics)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	inout		Initial conditions.
file_parameters	type(file_ini)	in		Simulation parameters ini file handler.
physics	type(flume_physics_object)	in		Physics (for the regions state conversion).

Call graph

set_equilibrium

Set the hydrostatic equilibrium of the well-balanced mode on the blocks interior (issue #73, M7-G2b): the atmosphere of the initial conditions at rest (zero velocity, region 1's field, no pulse or seed); the ghost cells are filled by the residual's ghost update. Fatal if the initial conditions are not an atmosphere.

subroutine set_equilibrium(self, field, q)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	in		Initial conditions.
field	type(field_object)	in		Field (realm component, threaded in).
q	real(kind=R8P)	inout		Conservative variables.

Call graph

set_initial_conditions

Set initial conditions on the blocks interior; ghost cells are filled by the following ghost update.

subroutine set_initial_conditions(self, field, q)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	in		Initial conditions.
field	type(field_object)	in		Field (realm component, threaded in).
q	real(kind=R8P)	inout		Conservative variables.

Call graph

hash_cell

Return two pseudo-random numbers in [-1, 1) that depend only on the block Morton code and the cell indexes.

Integer mixing with no integer overflow, so the same bits on every compiler; the cell key never depends on the rank that owns the block. The key is chained (mix the code, then fold in i, j, k one at a time and mix again), so keys whose bits overlap, such as small codes and small indexes, do not collide. Xorshift alone is linear over GF(2), which would make the numbers an affine function of the key bits (exactly equidistributed, spatially structured): each mix also multiplies by an odd constant modulo 2^52 (xorshift* idea, Vigna 2016), in 26-bit limbs to stay in range.

Attributes: pure

subroutine hash_cell(code, i, j, k, h)

Arguments

Name	Type	Intent	Attributes	Description
code	integer(kind=I8P)	in		Block Morton code.
i	integer(kind=I4P)	in		Cell indexes.
j	integer(kind=I4P)	in		Cell indexes.
k	integer(kind=I4P)	in		Cell indexes.
h	real(kind=R8P)	out		Pseudo-random numbers in [-1, 1).

Call graph

linear_wave_eigenvector

Return the right eigenvector of the right-going wave of the conservative MHD system at the primitive state prim (global frame), for the wave normal at angle degrees in the x-y plane: the state is rotated into the frame of the normal, R is the column of mhd_eigenvectors (x direction of that frame; waves ordered u_n - c_f, u_n - c_a, u_n - c_s, u_n, u_n + c_s, u_n + c_a, u_n + c_f, B_n), and its momentum and field are rotated back to x, y, z.

Attributes: pure

subroutine linear_wave_eigenvector(gamma, R, prim, wave, angle, r_global)

Arguments

Name	Type	Intent	Attributes	Description
gamma	real(kind=R8P)	in		Specific heats ratio.
R	real(kind=R8P)	in		Gas constant.
prim	real(kind=R8P)	in		Primitive state, global frame (r, u, v, w, p, bx, by, bz).
wave	character(len=*)	in		Wave family: fast, alfven, slow, entropy.
angle	real(kind=R8P)	in		Wave normal angle in the x-y plane (degrees).
r_global	real(kind=R8P)	out		Right eigenvector, global frame.

Call graph

linear_wave_eigenvector_inert

Return the right eigenvector of the right-going wave of a barotropic MHD model (issue #100, P5), as linear_wave_eigenvector: the column of the isothermal or barotropic mhd_eigenvectors (Kim et al. 1999; waves ordered u_n - c_f, u_n - c_a, u_n - c_s, u_n + c_s, u_n + c_a, u_n + c_f, then the inert slot and B_n), its momentum and field rotated back to x, y, z. There is no entropy mode (refused by the loader).

subroutine linear_wave_eigenvector_inert(model, eos_arg, prim, wave, angle, r_global)

Arguments

Name	Type	Intent	Attributes	Description
model	integer(kind=I4P)	in		Physical model id.
eos_arg	real(kind=R8P)	in		EOS parameter: c0^2, or the barotropic table id.
prim	real(kind=R8P)	in		Primitive state, global frame (r, u, v, w, p, bx, by, bz).
wave	character(len=*)	in		Wave family: fast, alfven, slow.
angle	real(kind=R8P)	in		Wave normal angle in the x-y plane (degrees).
r_global	real(kind=R8P)	out		Right eigenvector, global frame.

Call graph

mhd_vortex

Return the primitive state of the magnetised vortex at (x, y) (formulas in the module documentation).

Attributes: pure

subroutine mhd_vortex(prim0, mvortex, x, y, prim)

Arguments

Name	Type	Intent	Attributes	Description
prim0	real(kind=R8P)	in		Free stream primitive state (r, u, v, w, p, bx, by, bz), bx = by = 0.
mvortex	real(kind=R8P)	in		Vortex x0, y0, radius, kappa, mu.
x	real(kind=R8P)	in		Point coordinates.
y	real(kind=R8P)	in		Point coordinates.
prim	real(kind=R8P)	out		Primitive state.

Call graph

mhd_rotor

Return the primitive state of the MHD rotor at (x, y) (formulas in the module documentation).

Attributes: pure

subroutine mhd_rotor(prim0, rotor, x, y, prim)

Arguments

Name	Type	Intent	Attributes	Description
prim0	real(kind=R8P)	in		Ambient primitive state (r, u, v, w, p, bx, by, bz).
rotor	real(kind=R8P)	in		Rotor x0, y0, r0, r1, rho_in, v0.
x	real(kind=R8P)	in		Point coordinates.
y	real(kind=R8P)	in		Point coordinates.
prim	real(kind=R8P)	out		Primitive state.

Call graph

atmosphere_state

Primitive state of the atmosphere at rest at the point center (issue #73): isothermal, rho = rho_1 exp(g . x / T_1), p = rho T_1; polytropic, rho = rho_1 s^(1/(G-1)), p = p_1 s^(G/(G-1)), s = 1 + (G - 1) / G rho_1 / p_1 g . x; Rayleigh-Taylor, the two layers of the module header, with int_0^z dz / (a + b tanh z) = (a z - b ln(a cosh z + b sinh z) + b ln a) / (a^2 - b^2); zero velocity, region 1's field. Fatal where the polytrope has no mass.

subroutine atmosphere_state(self, center, prim)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	in		Initial conditions.
center	real(kind=R8P)	in		Point.
prim	real(kind=R8P)	out		Primitive state (first nprim used).

Call graph

rotate_to_global

Rotate the velocity and (MHD) the field of a primitive state from the frame of normal (normal, tangential, z) to x, y, z: v = v_n n + v_t t, n = normal / |normal|, t = (-n_y, n_x).

Attributes: pure

subroutine rotate_to_global(normal, nprim, prim)

Arguments

Name	Type	Intent	Attributes	Description
normal	real(kind=R8P)	in		Normal, (x, y), not normalised.
nprim	integer(kind=I4P)	in		Primitive keys number: 5 (Euler) or 8 (MHD).
prim	real(kind=R8P)	inout		Primitive state (r, u, v, w, p, bx, by, bz).

Call graph

set_isotropic_turbulence

Set the isotropic turbulence field (formulas in the module header): every block sums the modes on its own cells, the exponentials tabulated per axis.

subroutine set_isotropic_turbulence(self, field, q)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	in		The initial conditions.
field	type(field_object)	in		Field (cell centres).
q	real(kind=R8P)	inout		Conservative variables.

Call graph

isotropic_turbulence_modes

Return the modes of the isotropic turbulence field (formulas in the module header): the half space kx > 0, or kx = 0 and ky > 0, or kx = ky = 0 and kz > 0, swept in one order on every rank.

subroutine isotropic_turbulence_modes(hit, seed, kmax, mode_k, mode_a)

Arguments

Name	Type	Intent	Attributes	Description
hit	real(kind=R8P)	in		Density, u_rms, turbulent Mach number, k_0.
seed	integer(kind=I4P)	in		Seed of the phases.
kmax	integer(kind=I4P)	in		Largest wavenumber.
mode_k	integer(kind=I4P)	out	allocatable	Wavevectors, (3, mode).
mode_a	complex(kind=R8P)	out	allocatable	Complex amplitudes, (3, mode).

Call graph

channel_state

Return the primitive state (r, u, v, w, p) of a plane channel between y = +-h: uniform density and pressure, the laminar profile 3/2 U (1 - eta^2), and a perturbation of amplitude ch_pert U of a deterministic pattern, vanishing at the walls (2 waves along x, 6 along z), so identical on any rank count.

Attributes: pure

subroutine channel_state(gamma, ch, x, y, z, prim)

Arguments

Name	Type	Intent	Attributes	Description
gamma	real(kind=R8P)	in		Specific heats ratio.
ch	real(kind=R8P)	in		Density, bulk velocity, wall Mach number, half height, Lx, Lz, perturbation.
x	real(kind=R8P)	in		Cell centre.
y	real(kind=R8P)	in		Cell centre.
z	real(kind=R8P)	in		Cell centre.
prim	real(kind=R8P)	out		Primitive state.

Call graph

taylor_green

Return the primitive state (r, u, v, w, p) of the compressible Taylor-Green vortex at (x, y, z) (formulas in the module header).

Attributes: pure

subroutine taylor_green(gamma, tg, x, y, z, prim)

Arguments

Name	Type	Intent	Attributes	Description
gamma	real(kind=R8P)	in		Specific heats ratio.
tg	real(kind=R8P)	in		Density, velocity, Mach number, length.
x	real(kind=R8P)	in		Cell centre.
y	real(kind=R8P)	in		Cell centre.
z	real(kind=R8P)	in		Cell centre.
prim	real(kind=R8P)	out		Primitive state.

Call graph

isentropic_vortex

Return the conservative state of the isentropic vortex at (x, y) (formulas in the module documentation).

Attributes: pure

subroutine isentropic_vortex(gamma, prim, vortex, x, y, q)

Arguments

Name	Type	Intent	Attributes	Description
gamma	real(kind=R8P)	in		Specific heats ratio.
prim	real(kind=R8P)	in		Free stream primitive state (r, u, v, w, p).
vortex	real(kind=R8P)	in		Vortex x0, y0, radius, strength.
x	real(kind=R8P)	in		Point coordinates.
y	real(kind=R8P)	in		Point coordinates.
q	real(kind=R8P)	out		Conservative variables.

Call graph

becker_shock

Set Becker's exact viscous shock (Becker 1922, Z. Phys. 8, 321): steady 1-D Navier-Stokes along x, constant mu, Prandtl number 3/4 (the total enthalpy H = cp T + u^2 / 2 is then uniform), upstream state prim1 moving at u1 = mach a1 towards +x, shock centred at x0 (where u = (u1 + u2) / 2). With m = rho1 u1, the momentum integral is 4/3 mu u u' = m (gamma+1) / (2 gamma) (u - u1) (u - u2), whose solution is the implicit K (x - x0) = F(u) - F((u1 + u2) / 2), F(u) = A ln(u1 - u) + B ln(u - u2), K = 3 m (gamma+1) / (8 gamma mu), A = u1 / (u1 - u2), B = -u2 / (u1 - u2); F decreases monotonically, so bisection finds u.

subroutine becker_shock(gamma, R, cp, mu, prim1, mach, x0, x, q)

Arguments

Name	Type	Intent	Attributes	Description
gamma	real(kind=R8P)	in		Specific heats ratio, gas constant, cp.
R	real(kind=R8P)	in		Specific heats ratio, gas constant, cp.
cp	real(kind=R8P)	in		Specific heats ratio, gas constant, cp.
mu	real(kind=R8P)	in		Viscosity.
prim1	real(kind=R8P)	in		Upstream primitive state (rho, p used; velocity from the Mach number).
mach	real(kind=R8P)	in		Upstream Mach number.
x0	real(kind=R8P)	in		Shock position.
x	real(kind=R8P)	in		Abscissa.
q	real(kind=R8P)	out		Conservative variables.

Call graph

Functions

is_barotropic_safe

Return true if the IC type is defined for a barotropic gas (issue #100): its density, velocity and scalars are inputs or gamma-free formulas, and the pressure, which the barotropic gas ignores, enters nothing else. Refused: the isentropic vortex (rho from the adiabat, ^(1/(gamma-1))); the atmospheres and the polytrope (gamma-law equilibria); the dissipation ICs with their module; the MHD vortex (a pressure-balanced equilibrium) and the Orszag-Tang vortex (its density gamma^2 / (4 pi) presumes gamma = 5/3). The MHD ICs whose density, velocity and field are inputs or gamma-free formulas are accepted with the barotropic MHD models (P5); the linear wave takes their eigenvectors. The isothermal spheres (P4), singular and Bonnor-Ebert, are defined for the isothermal gas only.

Attributes: pure

Returns: logical

function is_barotropic_safe(self) result(safe)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	in		Initial conditions.

Call graph

description

Return a pretty-formatted object description.

Returns: character(len=:)

function description(self) result(desc)

Arguments

Name	Type	Intent	Attributes	Description
self	class(flume_ic_object)	in		Initial conditions.

Call graph

lane_emden_isothermal_table

Tabulate the isothermal Lane-Emden function (issue #100, P4, the Bonnor-Ebert sphere): psi'' + 2 psi' / xi = exp(-psi), psi(0) = psi'(0) = 0, at the nodes xi_n = n dxi up to xi_max and one beyond: the first node from the series psi = xi^2/6 - xi^4/120 + xi^6/1890 (the equation is singular at 0), the rest by classical RK4.

Attributes: pure

Returns: real(kind=R8P)

function lane_emden_isothermal_table(xi_max, dxi) result(tab)

Arguments

Name	Type	Intent	Attributes	Description
xi_max	real(kind=R8P)	in		Largest xi to cover.
dxi	real(kind=R8P)	in		Node spacing.

Call graph

lane_emden_isothermal_psi

Return psi(xi) from the Lane-Emden table by cubic Hermite interpolation of its nodes (values and slopes).

Attributes: pure

Returns: real(kind=R8P)

function lane_emden_isothermal_psi(tab, dxi, xi) result(psi)

Arguments

Name	Type	Intent	Attributes	Description
tab	real(kind=R8P)	in		psi and psi' at the nodes [2, 0:n].
dxi	real(kind=R8P)	in		Node spacing.
xi	real(kind=R8P)	in		Abscissa, within the table.

Call graph

inverse_square_mean

Mean of 1 / |x - c|^2 over the box [lo, hi] (issue #94, the singular isothermal sphere): by the divergence theorem, div((x - c) / |x - c|^2) = 1 / |x - c|^2 in 3D, the volume integral is the flux of (x - c) / |x - c|^2 through the six faces, int d / (d^2 + u^2 + v^2) du dv on a face at distance d from c: smooth unless the face passes through c, where it vanishes (d = 0). 8 x 8 Gauss-Legendre points per face.

Attributes: pure

Returns: real(kind=R8P)

function inverse_square_mean(lo, hi, c) result(mean)

Arguments

Name	Type	Intent	Attributes	Description
lo	real(kind=R8P)	in		Box corners.
hi	real(kind=R8P)	in		Box corners.
c	real(kind=R8P)	in		Singular point.

Call graph