# This file is part of COð˜•CEPT, the cosmological ð˜•-body code in Python.
# Copyright © 2015-2017 Jeppe Mosgaard Dakin.
#
# COð˜•CEPT is free software: You can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# COð˜•CEPT is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with COð˜•CEPT. If not, see http://www.gnu.org/licenses/
#
# The auther of COð˜•CEPT can be contacted at dakin(at)phys.au.dk
# The latest version of COð˜•CEPT is available at
# https://github.com/jmd-dk/concept/
# Import everything from the commons module.
# In the .pyx file, Cython declared variables will also get cimported.
from commons import *
# Entry point for the MacCormack method,
# which does time evolution of a fluid component.
@cython.header(# Arguments
component='Component',
ᔑdt='dict',
# Locals
attempt='int',
i='Py_ssize_t',
max_vacuum_corrections='int[::1]',
mc_step='int',
steps='Py_ssize_t[::1]',
Ï_ptr='double*',
Ïux_ptr='double*',
Ïuy_ptr='double*',
Ïuz_ptr='double*',
)
def maccormack(component, ᔑdt):
# Maximum allowed number of attempts to correct for
# negative densities, for the first and second MacCormack step.
max_vacuum_corrections = asarray([1, component.gridsize], dtype=C2np['int'])
# Extract fluid grid pointers
Ï_ptr = component.Ï.grid
Ïux_ptr = component.Ïux.grid
Ïuy_ptr = component.Ïuy.grid
Ïuz_ptr = component.Ïuz.grid
# Step/flux directions for the first MacCormack step
steps = next(maccormack_steps)
# The two MacCormack steps
for mc_step in range(2):
# Evolve the fluid variables. If this leads to negative
# densities, attempts are made to correct this.
for attempt in range(max_vacuum_corrections[mc_step]):
# Evolve fluid variables. In the first MacCormack step,
# the variables are re-evolved at each attempt. In the
# second MacCormack step, the variables should only be
# evolved once (vacuum correction may still take place
# multiple times).
if mc_step == 0 or attempt == 0:
# Nullify the starred grid buffers,
# so that they are ready to be populated
# by the following MacCormack step.
if mc_step == 0:
component.nullify_fluid_gridˣ()
# Compute starred variables from unstarred variables
evolve_fluid(component, ᔑdt, steps, mc_step)
# Nullify the Δ buffers, so that they are ready to
# be used by the following vacuum correction sweep.
component.nullify_fluid_Δ()
# Check and correct for density values heading dangerously
# fast towards negative values. If every density value
# is OK, accept this attempt at a MacCormack step as is.
if not correct_vacuum(component, mc_step):
break
else:
# None of the attempted MacCormack steps were accepted.
# If this is the second MacCormack step, this means that
# we have been unable to correct for negative densities.
if mc_step == 1:
abort('Giving up after {} failed attempts to remove negative densities in "{}"'
.format(max_vacuum_corrections[mc_step], component.name))
# Reverse step direction for the second MacCormack step
for i in range(3):
steps[i] *= -1
# The two MacCormack steps leave all values of all fluid variables
# with double their actual values. All grid values thus need
# to be halved. Note that no further communication is needed as we
# also halve the pseudo and ghost points.
for i in range(component.size):
Ï_ptr [i] *= 0.5
Ïux_ptr[i] *= 0.5
Ïuy_ptr[i] *= 0.5
Ïuz_ptr[i] *= 0.5
# Nullify the starred grid buffers and the Δ buffers,
# leaving these with no leftover junk.
component.nullify_fluid_gridˣ()
component.nullify_fluid_Δ()
# Infinite generator cycling through the 8 triples of
# step/flux directions, used in the maccormack function.
def generate_maccormack_steps():
steps = []
for sign in (+1, -1):
steps.append(sign*asarray((+1, +1, +1), dtype=C2np['Py_ssize_t']))
steps.append(sign*asarray((-1, +1, -1), dtype=C2np['Py_ssize_t']))
steps.append(sign*asarray((-1, -1, +1), dtype=C2np['Py_ssize_t']))
steps.append(sign*asarray((+1, -1, -1), dtype=C2np['Py_ssize_t']))
yield from itertools.cycle(steps)
maccormack_steps = generate_maccormack_steps()
# Function which evolve the fluid variables of a component
@cython.header(# Arguments
component='Component',
ᔑdt='dict',
steps='Py_ssize_t[::1]',
mc_step='int',
# Locals
h='double',
i='Py_ssize_t',
indices_local_start='Py_ssize_t[::1]',
indices_local_end='Py_ssize_t[::1]',
j='Py_ssize_t',
k='Py_ssize_t',
shape='tuple',
step_i='Py_ssize_t',
step_j='Py_ssize_t',
step_k='Py_ssize_t',
ux_ijk='double',
ux_sjk='double',
uy_ijk='double',
uy_isk='double',
uz_ijk='double',
uz_ijs='double',
w='double',
Ï='double[:, :, ::1]',
Ï_ijk='double',
Ï_ijs='double',
Ï_isk='double',
Ï_sjk='double',
Ïux='double[:, :, ::1]',
Ïux_ijk='double',
Ïux_ijs='double',
Ïux_isk='double',
Ïux_sjk='double',
ÏuxË£='double[:, :, ::1]',
Ïuy='double[:, :, ::1]',
Ïuy_ijk='double',
Ïuy_ijs='double',
Ïuy_isk='double',
Ïuy_sjk='double',
ÏuyË£='double[:, :, ::1]',
Ïuz='double[:, :, ::1]',
Ïuz_ijk='double',
Ïuz_ijs='double',
Ïuz_isk='double',
Ïuz_sjk='double',
ÏuzË£='double[:, :, ::1]',
ÏË£='double[:, :, ::1]',
)
def evolve_fluid(component, ᔑdt, steps, mc_step):
"""It is assumed that the unstarred and starred grids have
correctly populated pseudo and ghost points.
"""
# Physical grid spacing
h = boxsize/component.gridsize
# Exract steps in each direction
step_i, step_j, step_k = steps
# Arrays of start and end indices for the local part of the
# fluid grids, meaning disregarding pseudo points and ghost points.
# We have 2 ghost points in the beginning and 1 pseudo point and
# 2 ghost points in the end.
shape = component.shape
indices_local_start = asarray((2, 2, 2), dtype=C2np['Py_ssize_t'])
indices_local_end = asarray(shape , dtype=C2np['Py_ssize_t']) - 2 - 1
# Extract fluid grids
Ï = component.Ï .grid_mv
Ïux = component.Ïux.grid_mv
Ïuy = component.Ïuy.grid_mv
Ïuz = component.Ïuz.grid_mv
# Extract starred fluid grids
ÏË£ = component.Ï .gridË£_mv
ÏuxË£ = component.Ïux.gridË£_mv
ÏuyË£ = component.Ïuy.gridË£_mv
ÏuzË£ = component.Ïuz.gridË£_mv
# Get the equation of state parameter w at this instance in time
w = component.w()
# In the case of the second MacCormack step, the role of the
# starred and the unstarred variables should be swapped.
if mc_step == 1:
Ï , ÏË£ = ÏË£ , Ï
Ïux, ÏuxË£ = ÏuxË£, Ïux
Ïuy, ÏuyË£ = ÏuyË£, Ïuy
Ïuz, ÏuzË£ = ÏuzË£, Ïuz
# Loop which update the parsed starred variables
# from the parsed unstarred variables.
for i in range(ℤ[indices_local_start[0]], ℤ[indices_local_end[0]]):
for j in range(ℤ[indices_local_start[1]], ℤ[indices_local_end[1]]):
for k in range(ℤ[indices_local_start[2]], ℤ[indices_local_end[2]]):
# Density at this point
Ï_ijk = Ï[i, j, k]
# Momentum density components at this point
Ïux_ijk = Ïux[i, j, k]
Ïuy_ijk = Ïuy[i, j, k]
Ïuz_ijk = Ïuz[i, j, k]
# Velocity components at this point
ux_ijk = Ïux_ijk/Ï_ijk
uy_ijk = Ïuy_ijk/Ï_ijk
uz_ijk = Ïuz_ijk/Ï_ijk
# Density at forward (backward) points
Ï_sjk = Ï[i + step_i, j , k ]
Ï_isk = Ï[i , j + step_j, k ]
Ï_ijs = Ï[i , j , k + step_k]
# Momentum density components at forward (backward) points
Ïux_sjk = Ïux[i + step_i, j , k ]
Ïux_isk = Ïux[i , j + step_j, k ]
Ïux_ijs = Ïux[i , j , k + step_k]
Ïuy_sjk = Ïuy[i + step_i, j , k ]
Ïuy_isk = Ïuy[i , j + step_j, k ]
Ïuy_ijs = Ïuy[i , j , k + step_k]
Ïuz_sjk = Ïuz[i + step_i, j , k ]
Ïuz_isk = Ïuz[i , j + step_j, k ]
Ïuz_ijs = Ïuz[i , j , k + step_k]
# Velocity components at forward (backward) points
ux_sjk = Ïux_sjk/Ï_sjk
uy_isk = Ïuy_isk/Ï_isk
uz_ijs = Ïuz_ijs/Ï_ijs
# Flux of Ï (Ï*u)
Ï_flux = (+ step_i*(Ïux_sjk - Ïux_ijk)
+ step_j*(Ïuy_isk - Ïuy_ijk)
+ step_k*(Ïuz_ijs - Ïuz_ijk)
)
# Flux of Ïux (Ïux*u)
Ïux_flux = (+ step_i*(Ïux_sjk*ux_sjk - Ïux_ijk*ux_ijk)
+ step_j*(Ïux_isk*uy_isk - Ïux_ijk*uy_ijk)
+ step_k*(Ïux_ijs*uz_ijs - Ïux_ijk*uz_ijk)
# Pressure term
+ step_i*â„[w/(1 + w)]*(Ï_sjk - Ï_ijk)
)
# Flux of Ïuy (Ïuy*u)
Ïuy_flux = (+ step_i*(Ïuy_sjk*ux_sjk - Ïuy_ijk*ux_ijk)
+ step_j*(Ïuy_isk*uy_isk - Ïuy_ijk*uy_ijk)
+ step_k*(Ïuy_ijs*uz_ijs - Ïuy_ijk*uz_ijk)
# Pressure term
+ step_j*â„[w/(1 + w)]*(Ï_isk - Ï_ijk)
)
# Flux of Ïuz (Ïuz*u)
Ïuz_flux = (+ step_i*(Ïuz_sjk*ux_sjk - Ïuz_ijk*ux_ijk)
+ step_j*(Ïuz_isk*uy_isk - Ïuz_ijk*uy_ijk)
+ step_k*(Ïuz_ijs*uz_ijs - Ïuz_ijk*uz_ijk)
# Pressure term
+ step_k*â„[w/(1 + w)]*(Ï_ijs - Ï_ijk)
)
# Update Ï
ÏË£[i, j, k] += (# Initial value
+ Ï_ijk
# Flux
- â„[ᔑdt['aâ»Â²']/h]*Ï_flux
)
# Update Ïux
ÏuxË£[i, j, k] += (# Initial value
+ Ïux_ijk
# Flux
- â„[ᔑdt['aâ»Â²']/h]*Ïux_flux
)
# Update Ïuy
ÏuyË£[i, j, k] += (# Initial value
+ Ïuy_ijk
# Flux
- â„[ᔑdt['aâ»Â²']/h]*Ïuy_flux
)
# Update Ïuz
ÏuzË£[i, j, k] += (# Initial value
+ Ïuz_ijk
# Flux
- â„[ᔑdt['aâ»Â²']/h]*Ïuz_flux
)
# Populate the pseudo and ghost points with the updated values.
# Depedendent on whether we are doing the first of second
# MacCormack step (mc_step), the updated grids are really the
# starred grids (first MacCormack step) or the
# unstarred grids (second MacCormack step)
if mc_step == 0:
component.communicate_fluid_gridsˣ(mode='populate')
else: # mc_step == 1
component.communicate_fluid_grids(mode='populate')
# Function which checks for imminent vacuum in a fluid component
# and does one sweep of vacuum corrections.
@cython.header(# Arguments
component='Component',
mc_step='int',
# Locals
dist2='Py_ssize_t',
fac_smoothing='double',
fac_time='double',
i='Py_ssize_t',
indices_local_start='Py_ssize_t[::1]',
indices_local_end='Py_ssize_t[::1]',
j='Py_ssize_t',
k='Py_ssize_t',
m='Py_ssize_t',
mi='Py_ssize_t',
mj='Py_ssize_t',
mk='Py_ssize_t',
n='Py_ssize_t',
ni='Py_ssize_t',
nj='Py_ssize_t',
nk='Py_ssize_t',
shape='tuple',
timespan='double',
vacuum_imminent='bint',
ΔÏ='double[:, :, ::1]',
ΔÏ_ptr='double*',
ΔÏux='double[:, :, ::1]',
ΔÏux_ptr='double*',
ΔÏuy='double[:, :, ::1]',
ΔÏuy_ptr='double*',
ΔÏuz='double[:, :, ::1]',
ΔÏuz_ptr='double*',
Ï='double[:, :, ::1]',
Ï_correction='double',
Ï_ijk='double',
Ï_ptr='double*',
Ïux='double[:, :, ::1]',
Ïux_correction='double',
Ïux_ptr='double*',
ÏuxË£='double[:, :, ::1]',
Ïuy='double[:, :, ::1]',
Ïuy_correction='double',
Ïuy_ptr='double*',
ÏuyË£='double[:, :, ::1]',
Ïuz='double[:, :, ::1]',
Ïuz_correction='double',
Ïuz_ptr='double*',
ÏuzË£='double[:, :, ::1]',
ÏË£='double[:, :, ::1]',
ÏË£_ijk='double',
returns='bint',
)
def correct_vacuum(component, mc_step):
"""This function will detect and correct for imminent vacuum in a
fluid component. The vacuum detection is done differently depending
on the MacCormack step (the parsed mc_step). For the first
MacCormack step, vacuum is considered imminent if a density below
the vacuum density, Ï_vacuum, will be reached within timespan
similiar time steps. For the second MacCormack step, vacuum is
considered imminent if the density is below the vacuum density.
The vacuum correction is done by smoothing all fluid variables in
the 3x3x3 neighbouring cells souronding the vacuum cell.
The smoothing between each pair of cells, call them (Ïi, Ïj),
is given by
Ïi += fac_smoothing*(Ïj - Ïi)/r²,
Ïj += fac_smoothing*(Ïi - Ïj)/r²,
where r is the distance between the cells in grid units.
Whether or not any vacuum corrections were made is returned
as the return value.
Experimentally, it has been found that when
max_vacuum_corrections[0] == 1,
the following values give good, stable results:
timespan = 30
fac_smoothing = 1.5/(6/1 + 12/2 + 8/3)
"""
# In the case of the first MacCormack step, consider vacuum to be
# imminent if a cell will reach the vacuum density after this many
# similar time steps. Should be at least 1.
timespan = 30
# Amount of smoohing to apply when vacuum is detected.
# A numerator of 0 implies no smoothing.
# A numerator of 1 implies that in the most extreme case,
# a vacuum cell will be replaced with a weighted average of its
# 26 neighbour cells (all of the original cell will be distributed
# among these neighbors).
fac_smoothing = â„[1.5/(6/1 + 12/2 + 8/3)]
# Arrays of start and end indices for the local part of the
# fluid grids, meaning disregarding pseudo points and ghost points.
# We have 2 ghost points in the beginning and 1 pseudo point and
# 2 ghost points in the end.
shape = component.shape
indices_local_start = asarray([2, 2, 2], dtype=C2np['Py_ssize_t'])
indices_local_end = asarray(shape , dtype=C2np['Py_ssize_t']) - 2 - 1
# Extract memory views and pointers to the fluid variables
Ï = component.Ï .grid_mv
Ï_ptr = component.Ï .grid
ÏË£ = component.Ï .gridË£_mv
ÏË£_ptr = component.Ï .gridË£
Î”Ï = component.Ï .Δ_mv
ΔÏ_ptr = component.Ï .Δ
Ïux = component.Ïux.grid_mv
Ïux_ptr = component.Ïux.grid
ÏuxË£ = component.Ïux.gridË£_mv
ÏuxË£_ptr = component.Ïux.gridË£
ΔÏux = component.Ïux.Δ_mv
ΔÏux_ptr = component.Ïux.Δ
Ïuy = component.Ïuy.grid_mv
Ïuy_ptr = component.Ïuy.grid
ÏuyË£ = component.Ïuy.gridË£_mv
ÏuyË£_ptr = component.Ïuy.gridË£
ΔÏuy = component.Ïuy.Δ_mv
ΔÏuy_ptr = component.Ïuy.Δ
Ïuz = component.Ïuz.grid_mv
Ïuz_ptr = component.Ïuz.grid
ÏuzË£ = component.Ïuz.gridË£_mv
ÏuzË£_ptr = component.Ïuz.gridË£
ΔÏuz = component.Ïuz.Δ_mv
ΔÏuz_ptr = component.Ïuz.Δ
# In the case of the second MacCormack step, the role of the
# starred and the unstarred variables should be swapped.
if mc_step == 1:
Ï , ÏË£ = ÏË£ , Ï
Ï_ptr , ÏË£_ptr = ÏË£_ptr , Ï_ptr
Ïux , ÏuxË£ = ÏuxË£ , Ïux
Ïux_ptr, ÏuxË£_ptr = ÏuxË£_ptr, Ïux_ptr
Ïuy , ÏuyË£ = ÏuyË£ , Ïuy
Ïuy_ptr, ÏuyË£_ptr = ÏuyË£_ptr, Ïuy_ptr
Ïuz , ÏuzË£ = ÏuzË£ , Ïuz
Ïuz_ptr, ÏuzË£_ptr = ÏuzË£_ptr, Ïuz_ptr
# Loop over the local domain and check and compute
# corrections for imminent vacuum.
vacuum_imminent = False
for i in range(ℤ[indices_local_start[0]], ℤ[indices_local_end[0]]):
for j in range(ℤ[indices_local_start[1]], ℤ[indices_local_end[1]]):
for k in range(ℤ[indices_local_start[2]], ℤ[indices_local_end[2]]):
# Unstarred and starred density at this point
Ï_ijk = Ï [i, j, k]
ÏË£_ijk = ÏË£[i, j, k]
# Check for imminent vacuum.
# After the first MacCormack step, vacuum is considered
# to be imminent if a density below the vacuum density,
# Ï_vacuum, will be reached within timespan similiar
# time steps. That is, vacuum is imminent if
# Ï + timespan*dÏ < Ï_vacuum,
# where dÏ is the change in Ï from the first MacCormack
# step, given by dÏ = ½(ÏË£ - Ï), where the factor ½ is
# due to ÏË£ really holding double the change,
# ÏË£ = Ï + 2*dÏ. Put together, this means that vacuum
# is imminent if
# ÏË£ + Ï*(2/timespan - 1) < 2/timespan*Ï_vacuum.
# After the second MacCormack step, vacuum is considered
# to be imminent only if the density is lower than the
# vacuum density, Ï_vacuum. Because the starred
# variables hold double their actual values,
# this corresponds to
# ÏË£_ijk < 2*Ï_vacuum.
if ( (mc_step == 0 and Ï_ijk*â„[2/timespan - 1] + ÏË£_ijk < â„[2/timespan*Ï_vacuum])
or (mc_step == 1 and ÏË£_ijk < â„[2*Ï_vacuum])
):
vacuum_imminent = True
# The amount of smoothing to apply depends upon
# how far into the future densities below the vacuum
# density will be reached.
if mc_step == 0:
# The number of time steps before densities
# lower than the vacuum density is given by
# Ï + timesteps*dÏ == Ï_vacuum, dÏ = ½(ÏË£ - Ï).
# --> timesteps = 2*(Ï - Ï_vacuum)/(Ï - ÏË£).
fac_time = 0.5*(Ï_ijk - ÏË£_ijk)/(Ï_ijk - Ï_vacuum)
else: # mc_step == 1
# The density is already lower
# than the vaccuum density.
fac_time = 1
# Loop over all cell pairs (m, n) in the 3x3x3 block
# souronding the vacuum cell and apply smoothing.
for m in range(27):
# 3D indices of m'th cell
mi = i + relative_neighbour_indices_i[m]
mj = j + relative_neighbour_indices_j[m]
mk = k + relative_neighbour_indices_k[m]
for n in range(m + 1, 27):
# 3D indices of n'th cell
ni = i + relative_neighbour_indices_i[n]
nj = j + relative_neighbour_indices_j[n]
nk = k + relative_neighbour_indices_k[n]
# Distance squared between the two cells,
# in grid units (1 ≤ dist2 ≤ 12).
dist2 = (ni - mi)**2 + (nj - mj)**2 + (nk - mk)**2
# Compute vacuum corrections
Ï_correction = (Ï [ni, nj, nk] - â„[Ï [mi, mj, mk]])*â„[fac_smoothing*fac_time/dist2]
Ïux_correction = (Ïux[ni, nj, nk] - â„[Ïux[mi, mj, mk]])*â„[fac_smoothing*fac_time/dist2]
Ïuy_correction = (Ïuy[ni, nj, nk] - â„[Ïuy[mi, mj, mk]])*â„[fac_smoothing*fac_time/dist2]
Ïuz_correction = (Ïuz[ni, nj, nk] - â„[Ïuz[mi, mj, mk]])*â„[fac_smoothing*fac_time/dist2]
# Store vacuum corrections
Î”Ï [mi, mj, mk] += Ï_correction
ΔÏux[mi, mj, mk] += Ïux_correction
ΔÏuy[mi, mj, mk] += Ïuy_correction
ΔÏuz[mi, mj, mk] += Ïuz_correction
Î”Ï [ni, nj, nk] -= Ï_correction
ΔÏux[ni, nj, nk] -= Ïux_correction
ΔÏuy[ni, nj, nk] -= Ïuy_correction
ΔÏuz[ni, nj, nk] -= Ïuz_correction
# If vacuum is imminent on any process, consider it to be
# imminent on every process.
vacuum_imminent = allreduce(vacuum_imminent, op=MPI.LOR)
if vacuum_imminent:
# Communicate contributions to local vacuum corrections
# residing on other processes.
component.communicate_fluid_Δ(mode='add contributions')
# Local Δ buffers now store final values. Populate pseudo
# and ghost points of Δ buffers.
component.communicate_fluid_Δ(mode='populate')
# Apply vacuum corrections. Note that no further communication
# is needed as we also apply vacuum corrections to the
# pseudo and ghost points.
for i in range(component.size):
Ï_ptr [i] += ΔÏ_ptr [i]
Ïux_ptr[i] += ΔÏux_ptr[i]
Ïuy_ptr[i] += ΔÏuy_ptr[i]
Ïuz_ptr[i] += ΔÏuz_ptr[i]
# The return value should indicate whether or not
# vacuum corrections have been carried out.
return vacuum_imminent
# 1D memory views of relative indices to the 27 neighbours of a cell
# (itself included). These are thus effectively mappings from
# 1D indices to 3D indices.
cython.declare(relative_neighbour_indices_i='Py_ssize_t[::1]',
relative_neighbour_indices_j='Py_ssize_t[::1]',
relative_neighbour_indices_k='Py_ssize_t[::1]',
)
relative_neighbour_indices = asarray([(i, j, k) for i in range(-1, 2)
for j in range(-1, 2)
for k in range(-1, 2)], dtype=C2np['Py_ssize_t'])
relative_neighbour_indices_i = asarray(relative_neighbour_indices[:, 0]).copy()
relative_neighbour_indices_j = asarray(relative_neighbour_indices[:, 1]).copy()
relative_neighbour_indices_k = asarray(relative_neighbour_indices[:, 2]).copy()