May I know the most elegant way to integrate mass flux across a boundary to obtain the mass flow rate? I would like to track this value at some given time-steps, similar to how the forces/moments are tracked.
There is currently no direct means in PyFR of accomplishing this. However, it should not be difficult to add; the force integrator plugin contains all of the required mechanics for integrating quantities over boundaries. Thus with a few minor modifications it can be used to integrate mass.
Regards, Freddie.
Thank you. Does something like the following look like a step in the right direction?
from collections import defaultdict
import numpy as np
from pyfr.mpiutil import get_comm_rank_root, mpi
from pyfr.plugins.base import BaseSolnPlugin, SurfaceMixin, init_csv
class FluidFluxPlugin(SurfaceMixin, BaseSolnPlugin):
name = 'fluidflux'
systems = ['ac-euler', 'ac-navier-stokes', 'euler', 'navier-stokes']
formulations = ['dual', 'std']
dimensions = [2, 3]
def __init__(self, intg, cfgsect, suffix):
super().__init__(intg, cfgsect, suffix)
comm, rank, root = get_comm_rank_root()
# Output frequency
self.nsteps = self.cfg.getint(cfgsect, 'nsteps')
# Check if the system is incompressible
self._ac = intg.system.name.startswith('ac')
# Viscous correction
self._viscorr = self.cfg.get('solver', 'viscosity-correction', 'none')
# Constant variables
self._constants = self.cfg.items_as('constants', float)
# Underlying elements class
self.elementscls = intg.system.elementscls
# Boundary to integrate over
bc = f'bcon_{suffix}_p{intg.rallocs.prank}'
# Get the mesh and elements
mesh, elemap = intg.system.mesh, intg.system.ele_map
# See which ranks have the boundary
bcranks = comm.gather(bc in mesh, root=root)
# The root rank needs to open the output file
if rank == root:
if not any(bcranks):
raise RuntimeError(f'Boundary {suffix} does not exist')
# CSV header
header = ['t', 'massflux'][:2]
# Open
self.outf = init_csv(self.cfg, cfgsect, ','.join(header))
# Interpolation matrices and quadrature weights
self._m0 = m0 = {}
self._qwts = qwts = defaultdict(list)
# If we have the boundary then process the interface
if bc in mesh:
# Element indices and associated face normals
eidxs = defaultdict(list)
norms = defaultdict(list)
for etype, eidx, fidx, flags in mesh[bc].tolist():
eles = elemap[etype]
itype, proj, norm = eles.basis.faces[fidx]
ppts, pwts = self._surf_quad(itype, proj, flags='s')
nppts = len(ppts)
# Get phyical normals
pnorm = eles.pnorm_at(ppts, [norm]*nppts)[:, eidx]
eidxs[etype, fidx].append(eidx)
norms[etype, fidx].append(pnorm)
if (etype, fidx) not in m0:
m0[etype, fidx] = eles.basis.ubasis.nodal_basis_at(ppts)
qwts[etype, fidx] = pwts
self._eidxs = {k: np.array(v) for k, v in eidxs.items()}
self._norms = {k: np.array(v) for k, v in norms.items()}
def __call__(self, intg):
# Return if no output is due
if intg.nacptsteps % self.nsteps:
return
# MPI info
comm, rank, root = get_comm_rank_root()
# Solution matrices indexed by element type
solns = dict(zip(intg.system.ele_types, intg.soln))
nvars = self.nvars
# Flux
f = np.zeros(1)
for etype, fidx in self._m0:
# Get the interpolation operator
m0 = self._m0[etype, fidx]
nfpts, nupts = m0.shape
# Extract the relevant elements from the solution
uupts = solns[etype][..., self._eidxs[etype, fidx]]
# Interpolate to the face
ufpts = m0 @ uupts.reshape(nupts, -1)
ufpts = ufpts.reshape(nfpts, nvars, -1)
ufpts = ufpts.swapaxes(0, 1)
# Compute the velocity vector
rho = 1.0 if self._ac else self.elementscls.con_to_pri(ufpts, self.cfg)[0]
vs = np.array(self.elementscls.con_to_pri(ufpts, self.cfg)[1:-1])
# Get the quadrature weights and normal vectors
qwts = self._qwts[etype, fidx]
norms = self._norms[etype, fidx]
# Do the quadrature
f[0] += np.einsum('i,im,jim,dik->',qwts, rho, vs, norms)
# Reduce and output if we're the root rank
if rank != root:
comm.Reduce(f, None, op=mpi.SUM, root=root)
else:
comm.Reduce(mpi.IN_PLACE, f, op=mpi.SUM, root=root)
# Write
print(intg.tcurr, *f.ravel(), sep=',', file=self.outf)
# Flush to disk
self.outf.flush()
Yes, although you’ll always want to double check the einsum expression as they can be very tricky. It might be easier just to do rho*vs and pass that into einsum (which will dot with the normal and do the quadrature).
Regards, Freddie.
Thanks. There was some issue with the previous code but the following should work:
from collections import defaultdict
import numpy as np
from pyfr.mpiutil import get_comm_rank_root, mpi
from pyfr.plugins.base import BaseSolnPlugin, SurfaceMixin, init_csv
class FluidFluxPlugin(SurfaceMixin, BaseSolnPlugin):
name = 'fluidflux'
systems = ['ac-euler', 'ac-navier-stokes', 'euler', 'navier-stokes']
formulations = ['dual', 'std']
dimensions = [2, 3]
def __init__(self, intg, cfgsect, suffix):
super().__init__(intg, cfgsect, suffix)
comm, rank, root = get_comm_rank_root()
# Output frequency
self.nsteps = self.cfg.getint(cfgsect, 'nsteps')
# Check if the system is incompressible
self._ac = intg.system.name.startswith('ac')
# Viscous correction
self._viscorr = self.cfg.get('solver', 'viscosity-correction', 'none')
# Constant variables
self._constants = self.cfg.items_as('constants', float)
# Underlying elements class
self.elementscls = intg.system.elementscls
# Boundary to integrate over
bc = f'bcon_{suffix}_p{intg.rallocs.prank}'
# Get the mesh and elements
mesh, elemap = intg.system.mesh, intg.system.ele_map
# See which ranks have the boundary
bcranks = comm.gather(bc in mesh, root=root)
# The root rank needs to open the output file
if rank == root:
if not any(bcranks):
raise RuntimeError(f'Boundary {suffix} does not exist')
# CSV header
header = ['t', 'massflux'][:2]
# Open
self.outf = init_csv(self.cfg, cfgsect, ','.join(header))
# Interpolation matrices and quadrature weights
self._m0 = m0 = {}
self._qwts = qwts = defaultdict(list)
# If we have the boundary then process the interface
if bc in mesh:
# Element indices and associated face normals
eidxs = defaultdict(list)
norms = defaultdict(list)
for etype, eidx, fidx, flags in mesh[bc].tolist():
eles = elemap[etype]
itype, proj, norm = eles.basis.faces[fidx]
ppts, pwts = self._surf_quad(itype, proj, flags='s')
nppts = len(ppts)
# Get phyical normals
pnorm = eles.pnorm_at(ppts, [norm]*nppts)[:, eidx]
eidxs[etype, fidx].append(eidx)
norms[etype, fidx].append(pnorm)
if (etype, fidx) not in m0:
m0[etype, fidx] = eles.basis.ubasis.nodal_basis_at(ppts)
qwts[etype, fidx] = pwts
self._eidxs = {k: np.array(v) for k, v in eidxs.items()}
self._norms = {k: np.array(v) for k, v in norms.items()}
def __call__(self, intg):
# Return if no output is due
if intg.nacptsteps % self.nsteps:
return
# MPI info
comm, rank, root = get_comm_rank_root()
# Solution matrices indexed by element type
solns = dict(zip(intg.system.ele_types, intg.soln))
nvars = self.nvars
# Flux
f = np.zeros(1)
for etype, fidx in self._m0:
# Get the interpolation operator
m0 = self._m0[etype, fidx]
nfpts, nupts = m0.shape
# Extract the relevant elements from the solution
uupts = solns[etype][..., self._eidxs[etype, fidx]]
# Interpolate to the face
ufpts = m0 @ uupts.reshape(nupts, -1)
ufpts = ufpts.reshape(nfpts, nvars, -1)
ufpts = ufpts.swapaxes(0, 1)
# Convert conservative variables to primitive variables
pri_vars = self.elementscls.con_to_pri(ufpts, self.cfg)
vs = np.array(pri_vars[1:-1]) # Velocity components
if self._ac:
# Artificial compressibility: density is not a primitive variable
# Set rho to 1 with the appropriate shape
rho = np.ones_like(vs[0])
else:
# Compressible flow: extract rho from primitive variables
rho = pri_vars[0]
# Compute rhovs: multiply rho and vs
rhovs = rho[None, :, :] * vs
# Get the quadrature weights and normal vectors
qwts = self._qwts[etype, fidx]
norms = self._norms[etype, fidx]
# Do the quadrature
f[0] += np.einsum('i,jim,mij->', qwts, rhovs, norms)
# Reduce and output if we're the root rank
if rank != root:
comm.Reduce(f, None, op=mpi.SUM, root=root)
else:
comm.Reduce(mpi.IN_PLACE, f, op=mpi.SUM, root=root)
# Write
print(intg.tcurr, *f.ravel(), sep=',', file=self.outf)
# Flush to disk
self.outf.flush()
This seems to work on my simple cases but I am facing some indirectly related issues (like for instance I am not able to maintain a mass flow rate through a boundary with an imposed velocity, perhaps due to weakly-imposed boundary conditions?). I will keep testing it further and do a pull request if I get the chance.
It is not always possible to maintain a mass flow rate through a boundary. Consider a channel with a high-subsonic inflow. Here, you should have no trouble getting the expected mass-flow rate at the boundary. Now, consider putting a body (such as a aerofoil) in the channel. As the flow approaches the aerofoil it will choke which, in turn, will limit the amount of mass which can enter the domain upstream. The mass flow rate will therefore be less than you expect.
Regards, Freddie.