The density of a condensed mixture¶

The density of a condensed mixture from its components' densities.

An ideal mixture's volumes add: a kilogram of it takes the mass-weighted sum of its components' volumes, 1 / rho = sum w_i / rho_i, w_i = m_i / sum m the mass fractions (no excess volume of mixing). Over a game's whole table of parts: each part's mixture, its mass (a fraction, or any masses: only their proportions count) and its component, a relation to whatever table the game's column relates to, whose density it reads (lusine.materials's mechanical_phases, or the game's own). It folds the parts of the item's own mixture; each sum is added in the fixed tree. A gas's density is lusine.properties.gas_density's. Source: Perry's Chemical Engineers' Handbook, 9th ed., section 2 (densities of mixtures, the additive-volume rule). Checked against a reference in numpy, which has its own tests.

Source

Perry's Chemical Engineers' Handbook, 9th ed., section 2 (densities of mixtures, the additive-volume rule).

Proof

reference, limit, invariant.

mixture_density.trait

namespace lusine.properties;

float<kg/m3> mixture_density(list<struct { relation mixture; float fraction; relation component; }> parts)
{
  let total = sum(fraction) where mixture == it;
  return total / (sum(fraction / component.density) where mixture == it);
}

An ideal mixture's volumes add, 1 / ρ = Σ w_i / ρ_i. The kernel takes a game's whole table of parts, list<struct { relation mixture; float fraction; relation component; }>: each part's mixture, its mass and its component, a relation to whatever table holds the component's density.

In [1]:

import os, sys
root = os.getcwd()
while not os.path.isfile(os.path.join(root, 'docsite.json')):
    root = os.path.dirname(root)
sys.path.insert(0, os.path.join(root, 'tools'))
sys.path.insert(0, os.path.join(root, 'kernels'))
import numpy as np
import matplotlib.pyplot as plt
from lusine import trait_run
F32 = np.float32
U = float(np.finfo(F32).eps) / 2        # binary32's unit roundoff, 2^-24

import os, sys root = os.getcwd() while not os.path.isfile(os.path.join(root, 'docsite.json')): root = os.path.dirname(root) sys.path.insert(0, os.path.join(root, 'tools')) sys.path.insert(0, os.path.join(root, 'kernels')) import numpy as np import matplotlib.pyplot as plt from lusine import trait_run F32 = np.float32 U = float(np.finfo(F32).eps) / 2 # binary32's unit roundoff, 2^-24

The kernel¶

The law as the library states it.

In [2]:

print(trait_run.source('properties/mixture_density.trait'))

print(trait_run.source('properties/mixture_density.trait'))

// lusine.properties.mixture_density: the density of a condensed mixture from its components' densities. An ideal
// mixture's volumes add: a kilogram of it takes the mass-weighted sum of its components' volumes, 1 / rho = sum w_i /
// rho_i, w_i = m_i / sum m the mass fractions (no excess volume of mixing). Over a game's whole table of parts: each
// part's mixture, its mass (a fraction, or any masses: only their proportions count) and its
// component, a relation to whatever table the game's column relates to, whose density it reads (lusine.materials's
// mechanical_phases, or the game's own). It folds the parts of the item's own mixture; each sum is added in the fixed
// tree. A gas's density is lusine.properties.gas_density's.
// Source: Perry's Chemical Engineers' Handbook, 9th ed., section 2 (densities of mixtures, the additive-volume rule).
// Checked against a reference in numpy, which has its own tests.

namespace lusine.properties;

float<kg/m3> mixture_density(list<struct { relation mixture; float fraction; relation component; }> parts)
{
  let total = sum(fraction) where mixture == it;
  return total / (sum(fraction / component.density) where mixture == it);
}

The program¶

A game's tables, written out: its components (water, 1000 kg/m³, and quartz, 2650), its mixtures and its parts, each a mixture's mass of a component: equal masses of the two (wet sand), water alone in 5 kg, three parts water to one of quartz in kilograms and in hundreds of kilograms, and quartz's mass fraction from 0 to 1 in twentieths. main() binds the kernel over the mixtures; the planted fault below binds the page's own wrong kernel in its place.

In [3]:

DENSITY = {'water': 1000.0, 'quartz': 2650.0}
SWEEP = [i / 20 for i in range(21)]
MIXTURES = {'wet': {'water': 1.0, 'quartz': 1.0}, 'water': {'water': 5.0}, 'three': {'water': 3.0, 'quartz': 1.0},
            'three00': {'water': 300.0, 'quartz': 100.0}}
MIXTURES.update({f'sweep{i}': {c: m for c, m in (('water', 1 - w), ('quartz', w)) if m > 0}
                 for i, w in enumerate(SWEEP)})


def program(kernel='lusine.properties.mixture_density', extra=''):
    components = ', '.join(f'{{ name = "{c}", density = {d!r} kg/m3 }}' for c, d in DENSITY.items())
    mixtures = ', '.join(f'{{ name = "{m}" }}' for m in MIXTURES)
    parts = ',\n'.join(f'  {{ mixture = mixtures["{m}"], fraction = {float(x)!r}, component = components["{c}"] }}'
                       for m, cs in MIXTURES.items() for c, x in cs.items())
    return f'''namespace sand;

lusine.version version = {{ {{ number = 1 }} }};

struct component {{ text name; float<kg/m3> density; }};
table<component> components(capacity = 2, key = name) = {{ {components} }};
struct mixture {{ text name; float<kg/m3> density; }};
table<mixture> mixtures(capacity = 32, key = name) = {{ {mixtures} }};
struct part {{ relation<mixtures> mixture; float fraction; relation<components> component; }};
table<part> parts(capacity = 64) = {{
{parts},
}};
{extra}
void main() {{ mixtures.density = {kernel}(parts); }}
'''

LAW = program()
print(LAW)


def densities(text=LAW):
    run = trait_run.Program(['properties/mixture_density.trait'], text, {})
    return dict(zip(MIXTURES, run.steps('sand.main', 1, [('sand.mixtures', 'density')])[0][0]))

R = densities()

DENSITY = {'water': 1000.0, 'quartz': 2650.0} SWEEP = [i / 20 for i in range(21)] MIXTURES = {'wet': {'water': 1.0, 'quartz': 1.0}, 'water': {'water': 5.0}, 'three': {'water': 3.0, 'quartz': 1.0}, 'three00': {'water': 300.0, 'quartz': 100.0}} MIXTURES.update({f'sweep{i}': {c: m for c, m in (('water', 1 - w), ('quartz', w)) if m > 0} for i, w in enumerate(SWEEP)}) def program(kernel='lusine.properties.mixture_density', extra=''): components = ', '.join(f'{{ name = "{c}", density = {d!r} kg/m3 }}' for c, d in DENSITY.items()) mixtures = ', '.join(f'{{ name = "{m}" }}' for m in MIXTURES) parts = ',\n'.join(f' {{ mixture = mixtures["{m}"], fraction = {float(x)!r}, component = components["{c}"] }}' for m, cs in MIXTURES.items() for c, x in cs.items()) return f'''namespace sand; lusine.version version = {{ {{ number = 1 }} }}; struct component {{ text name; float<kg/m3> density; }}; table<component> components(capacity = 2, key = name) = {{ {components} }}; struct mixture {{ text name; float<kg/m3> density; }}; table<mixture> mixtures(capacity = 32, key = name) = {{ {mixtures} }}; struct part {{ relation<mixtures> mixture; float fraction; relation<components> component; }}; table<part> parts(capacity = 64) = {{ {parts}, }}; {extra} void main() {{ mixtures.density = {kernel}(parts); }} ''' LAW = program() print(LAW) def densities(text=LAW): run = trait_run.Program(['properties/mixture_density.trait'], text, {}) return dict(zip(MIXTURES, run.steps('sand.main', 1, [('sand.mixtures', 'density')])[0][0])) R = densities()

namespace sand;

lusine.version version = { { number = 1 } };

struct component { text name; float<kg/m3> density; };
table<component> components(capacity = 2, key = name) = { { name = "water", density = 1000.0 kg/m3 }, { name = "quartz", density = 2650.0 kg/m3 } };
struct mixture { text name; float<kg/m3> density; };
table<mixture> mixtures(capacity = 32, key = name) = { { name = "wet" }, { name = "water" }, { name = "three" }, { name = "three00" }, { name = "sweep0" }, { name = "sweep1" }, { name = "sweep2" }, { name = "sweep3" }, { name = "sweep4" }, { name = "sweep5" }, { name = "sweep6" }, { name = "sweep7" }, { name = "sweep8" }, { name = "sweep9" }, { name = "sweep10" }, { name = "sweep11" }, { name = "sweep12" }, { name = "sweep13" }, { name = "sweep14" }, { name = "sweep15" }, { name = "sweep16" }, { name = "sweep17" }, { name = "sweep18" }, { name = "sweep19" }, { name = "sweep20" } };
struct part { relation<mixtures> mixture; float fraction; relation<components> component; };
table<part> parts(capacity = 64) = {
  { mixture = mixtures["wet"], fraction = 1.0, component = components["water"] },
  { mixture = mixtures["wet"], fraction = 1.0, component = components["quartz"] },
  { mixture = mixtures["water"], fraction = 5.0, component = components["water"] },
  { mixture = mixtures["three"], fraction = 3.0, component = components["water"] },
  { mixture = mixtures["three"], fraction = 1.0, component = components["quartz"] },
  { mixture = mixtures["three00"], fraction = 300.0, component = components["water"] },
  { mixture = mixtures["three00"], fraction = 100.0, component = components["quartz"] },
  { mixture = mixtures["sweep0"], fraction = 1.0, component = components["water"] },
  { mixture = mixtures["sweep1"], fraction = 0.95, component = components["water"] },
  { mixture = mixtures["sweep1"], fraction = 0.05, component = components["quartz"] },
  { mixture = mixtures["sweep2"], fraction = 0.9, component = components["water"] },
  { mixture = mixtures["sweep2"], fraction = 0.1, component = components["quartz"] },
  { mixture = mixtures["sweep3"], fraction = 0.85, component = components["water"] },
  { mixture = mixtures["sweep3"], fraction = 0.15, component = components["quartz"] },
  { mixture = mixtures["sweep4"], fraction = 0.8, component = components["water"] },
  { mixture = mixtures["sweep4"], fraction = 0.2, component = components["quartz"] },
  { mixture = mixtures["sweep5"], fraction = 0.75, component = components["water"] },
  { mixture = mixtures["sweep5"], fraction = 0.25, component = components["quartz"] },
  { mixture = mixtures["sweep6"], fraction = 0.7, component = components["water"] },
  { mixture = mixtures["sweep6"], fraction = 0.3, component = components["quartz"] },
  { mixture = mixtures["sweep7"], fraction = 0.65, component = components["water"] },
  { mixture = mixtures["sweep7"], fraction = 0.35, component = components["quartz"] },
  { mixture = mixtures["sweep8"], fraction = 0.6, component = components["water"] },
  { mixture = mixtures["sweep8"], fraction = 0.4, component = components["quartz"] },
  { mixture = mixtures["sweep9"], fraction = 0.55, component = components["water"] },
  { mixture = mixtures["sweep9"], fraction = 0.45, component = components["quartz"] },
  { mixture = mixtures["sweep10"], fraction = 0.5, component = components["water"] },
  { mixture = mixtures["sweep10"], fraction = 0.5, component = components["quartz"] },
  { mixture = mixtures["sweep11"], fraction = 0.44999999999999996, component = components["water"] },
  { mixture = mixtures["sweep11"], fraction = 0.55, component = components["quartz"] },
  { mixture = mixtures["sweep12"], fraction = 0.4, component = components["water"] },
  { mixture = mixtures["sweep12"], fraction = 0.6, component = components["quartz"] },
  { mixture = mixtures["sweep13"], fraction = 0.35, component = components["water"] },
  { mixture = mixtures["sweep13"], fraction = 0.65, component = components["quartz"] },
  { mixture = mixtures["sweep14"], fraction = 0.30000000000000004, component = components["water"] },
  { mixture = mixtures["sweep14"], fraction = 0.7, component = components["quartz"] },
  { mixture = mixtures["sweep15"], fraction = 0.25, component = components["water"] },
  { mixture = mixtures["sweep15"], fraction = 0.75, component = components["quartz"] },
  { mixture = mixtures["sweep16"], fraction = 0.19999999999999996, component = components["water"] },
  { mixture = mixtures["sweep16"], fraction = 0.8, component = components["quartz"] },
  { mixture = mixtures["sweep17"], fraction = 0.15000000000000002, component = components["water"] },
  { mixture = mixtures["sweep17"], fraction = 0.85, component = components["quartz"] },
  { mixture = mixtures["sweep18"], fraction = 0.09999999999999998, component = components["water"] },
  { mixture = mixtures["sweep18"], fraction = 0.9, component = components["quartz"] },
  { mixture = mixtures["sweep19"], fraction = 0.050000000000000044, component = components["water"] },
  { mixture = mixtures["sweep19"], fraction = 0.95, component = components["quartz"] },
  { mixture = mixtures["sweep20"], fraction = 1.0, component = components["quartz"] },
};

void main() { mixtures.density = lusine.properties.mixture_density(parts); }

Reference: wet sand¶

Equal masses of water and quartz: their volumes added, 1452.05 kg/m³.

In [4]:

print(f"{R['wet']:.3f} kg/m3")
assert abs(R['wet'] - 1452.05) < 0.01

print(f"{R['wet']:.3f} kg/m3") assert abs(R['wet'] - 1452.05) < 0.01

1452.055 kg/m3

Reference: bit for bit, the folds in binary32¶

Every mixture is a binary32 model's bit for bit: its parts in slot order, their masses and their mass / ρ_i each added in Trait's fixed tree, and the first divided by the second. The float64 reference agrees within 4 × 2⁻²⁴.

In [5]:

from lusine.kernels.properties.mixture_density import mixture_density as reference


def tree(values):
    if not values:
        return F32(0.0)
    if len(values) == 1:
        return F32(values[0])
    mid = (len(values) + 1) // 2
    return F32(tree(values[:mid]) + tree(values[mid:]))


def model(composition):
    total = tree([F32(m) for m in composition.values()])
    return F32(total / tree([F32(F32(m) / F32(DENSITY[c])) for c, m in composition.items()]))


def model_holds(R):
    return all(F32(R[m]) == model(c) for m, c in MIXTURES.items())

assert model_holds(R)
worst = max(abs(R[m] - reference(c, DENSITY)) / R[m] for m, c in MIXTURES.items())
assert worst <= 4 * U, worst / U
print(f'{len(MIXTURES)} mixtures bit for bit with the binary32 model; the float64 reference within {worst / U:.1f} x 2^-24')

from lusine.kernels.properties.mixture_density import mixture_density as reference def tree(values): if not values: return F32(0.0) if len(values) == 1: return F32(values[0]) mid = (len(values) + 1) // 2 return F32(tree(values[:mid]) + tree(values[mid:])) def model(composition): total = tree([F32(m) for m in composition.values()]) return F32(total / tree([F32(F32(m) / F32(DENSITY[c])) for c, m in composition.items()])) def model_holds(R): return all(F32(R[m]) == model(c) for m, c in MIXTURES.items()) assert model_holds(R) worst = max(abs(R[m] - reference(c, DENSITY)) / R[m] for m, c in MIXTURES.items()) assert worst <= 4 * U, worst / U print(f'{len(MIXTURES)} mixtures bit for bit with the binary32 model; the float64 reference within {worst / U:.1f} x 2^-24')

25 mixtures bit for bit with the binary32 model; the float64 reference within 1.2 x 2^-24

Limit: one component alone¶

A mixture of one component is that component's density, within an ulp (1 / fl(1 / ρ)), at any mass.

In [6]:

for m in ('water', 'sweep0', 'sweep20'):
    want = DENSITY['quartz' if m == 'sweep20' else 'water']
    assert abs(R[m] - want) <= 2 * U * want, (m, R[m])
print('water alone', R['water'], '; quartz alone', R['sweep20'])

for m in ('water', 'sweep0', 'sweep20'): want = DENSITY['quartz' if m == 'sweep20' else 'water'] assert abs(R[m] - want) <= 2 * U * want, (m, R[m]) print('water alone', R['water'], '; quartz alone', R['sweep20'])

water alone 1000.0 ; quartz alone 2650.0

Invariant: proportions, between the components, and the volumes add¶

The same composition in kilograms or hundreds of kilograms has the same density; across every mass fraction the mixture lies between its components, rising with the quartz; and a kilogram's volume is the components' volumes added, within the sums' rounding.

In [7]:

def volumes_hold(R):
    rho = np.array([R[f'sweep{i}'] for i in range(len(SWEEP))])
    w = np.array(SWEEP)
    return (bool(np.all(np.abs(1 / rho - ((1 - w) / 1000.0 + w / 2650.0)) <= 6 * U / rho))
            and abs(R['three'] - R['three00']) <= 4 * U * R['three'])

rho = np.array([R[f'sweep{i}'] for i in range(len(SWEEP))])
assert volumes_hold(R) and np.all((rho >= F32(1000.0) * (1 - 2 * U)) & (rho <= 2650.0 * (1 + 2 * U)))
assert np.all(np.diff(rho) > 0)
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(SWEEP, rho, 'o-')
ax.set_xlabel('quartz, mass fraction')
ax.set_ylabel('density (kg/m3)')
plt.show()

def volumes_hold(R): rho = np.array([R[f'sweep{i}'] for i in range(len(SWEEP))]) w = np.array(SWEEP) return (bool(np.all(np.abs(1 / rho - ((1 - w) / 1000.0 + w / 2650.0)) <= 6 * U / rho)) and abs(R['three'] - R['three00']) <= 4 * U * R['three']) rho = np.array([R[f'sweep{i}'] for i in range(len(SWEEP))]) assert volumes_hold(R) and np.all((rho >= F32(1000.0) * (1 - 2 * U)) & (rho <= 2650.0 * (1 + 2 * U))) assert np.all(np.diff(rho) > 0) fig, ax = plt.subplots(figsize=(7, 3)) ax.plot(SWEEP, rho, 'o-') ax.set_xlabel('quartz, mass fraction') ax.set_ylabel('density (kg/m3)') plt.show()

Planted fault¶

A wrong kernel over the same list, written here in Trait and bound in the library's place: the densities' mass-weighted mean, as if the densities added rather than the volumes. Wet sand, the volumes and the binary32 model each catch it.

In [8]:

MEAN = program('mean', '''
float<kg/m3> mean(list<struct { relation mixture; float fraction; relation component; }> parts)
{
  let total = sum(fraction) where mixture == it;
  return (sum(fraction * component.density) where mixture == it) / total;
}
''')
W = densities(MEAN)
assert abs(W['wet'] - 1452.05) > 0.01 and not volumes_hold(W) and not model_holds(W)
print(f"PLANTED: the densities' mean gives wet sand {W['wet']:.1f} kg/m3, and is caught")

MEAN = program('mean', ''' float<kg/m3> mean(list<struct { relation mixture; float fraction; relation component; }> parts) { let total = sum(fraction) where mixture == it; return (sum(fraction * component.density) where mixture == it) / total; } ''') W = densities(MEAN) assert abs(W['wet'] - 1452.05) > 0.01 and not volumes_hold(W) and not model_holds(W) print(f"PLANTED: the densities' mean gives wet sand {W['wet']:.1f} kg/m3, and is caught")

PLANTED: the densities' mean gives wet sand 1825.0 kg/m3, and is caught