Extending Mocca

Most studies need only the Reference functions. This page is for adding your own kinds of step and unit, and for working with the simulation model that a Process builds. That model is a Jutul MultiModel, with one model per bed and device, and the steps of the cycle become its forces.

Adding a kind of step

The set of steps is open: define a subtype of AbstractStep and a method of step_settings.

Mocca.Interface.AbstractStep Type

AbstractStep

A step of a cycle, written in process terms. The built-in steps are Pressurise, Adsorb, Blowdown, Evacuate, Purge, Equalise and Idle.

To add a kind of step, define a subtype with a duration field (and usually bed and name fields) and a method of step_settings. Optionally add methods of end_pressure (if the step leaves its bed at a known pressure), step_beds (if it acts on several beds) and default_step_name. Register it with register_step! to read and write it in files; files give its fields as keywords. For example, a rinse of the bed with heavy product at the feed pressure:

Base.@kwdef mutable struct Rinse <: Mocca.AbstractStep
    duration::Float64
    bed::Union{Nothing, Symbol} = nothing
    name::String = "rinse"
end

function Mocca.step_settings(s::Rinse, ctx::Mocca.StepContext)
    d = Mocca.find_device(ctx, (:feed,))
    return [d => (law = Mocca.valve(ctx, d), inlet = Mocca.boundary(ctx, d, ctx.fs.feed.pressure))]
end
Mocca.end_pressure(s::Rinse, ctx) = ctx.fs.feed.pressure
Mocca.register_step!("rinse", Rinse)

source

Mocca.Interface.StepContext Type

StepContext

What a step sees when its settings are built: the flowsheet fs, its bed (resolved: a step that names no bed acts on the flowsheet's only bed), the bed's pressure at the start of the step (or nothing if no earlier step set it), the start pressures of all beds (pressures), the default ramp_rate [1/s] and a location for error messages. Helpers for step authors: find_device, valve, boundary, start_pressure and ramp.

source

Mocca.Interface.compile_cycle Function

compile_cycle(cycle, fs; ramp_rate = 0.5) → CompiledCycle

The stages of cycle on the flowsheet fs. Pressure ramps run at ramp_rate [1/s] unless a step sets its own.

source

Mocca.Interface.step_settings Function

step_settings(step, ctx::StepContext) → Vector{Pair{Symbol, NamedTuple}}

The devices step opens and their settings, as device => (law = ..., inlet = ..., outlet = ...); see Stage. Devices not listed are closed.

source

Mocca.Interface.end_pressure Function

end_pressure(step, ctx) → pressure or nothing

The pressure [Pa] step leaves its bed at, for the ramp of the next step, or nothing if it leaves it as it was.

source

Mocca.Interface.step_beds Function

step_beds(step) → Vector{Symbol}

The beds step acts on; by default its bed field, if set.

source

Mocca.Interface.default_step_name Function

Name of a step of this kind, when it is not given one.

source

Mocca.Interface.register_step! Function

register_step!(name, T)

Make the step type T readable from and writable to process files under "type": name. Its fields are given as keywords.

source

Mocca.Interface.find_device Function

find_device(ctx, roles; port = nothing, device = nothing) → Symbol

The device of the step's bed with the first of roles that is present, attached at port (:bottom or :top) if given. With device, check that it is attached to the bed and return it.

source

Mocca.Interface.valve Function

valve(ctx, device; conductance = nothing) → LinearValve

A fully open valve for device, with its default_conductance unless conductance [m³/(s·Pa)] is given.

source

Mocca.Interface.boundary Function

boundary(ctx, device, pressure) → PortCondition

The boundary of device with pressure (a number or a profile such as ExponentialRamp), and the temperature and composition of its default boundary.

source

Mocca.Interface.start_pressure Function

start_pressure(ctx, from) → Float64

from if given, else the pressure the bed was left at by the step before.

source

Mocca.Interface.ramp Function

ExponentialRamp from p0 to p1 at rate [1/s], by default the cycle's.

source

The layer underneath

The objects of the Reference build these, which can also be used directly. A Flowsheet holds a UnitNetwork: the unit models and the connections between their ports, from which the simulation model is built.

Mocca.Interface.fixed_bed_model Function

fixed_bed_model(gas, feed, adsorbent, column; ncells = 30, thermal = WithWall(),
    dispersion_velocity = 1.0, dispersion = nothing, permeability = nothing)

The model of a fixed bed of adsorbent packed in column, for gas, on ncells cells. The permeability (Kozeny–Carman, from the bed porosity and particle size) and the axial dispersion (0.7 D_m + 0.5 v d_p, at the interstitial dispersion_velocity [m/s]) are computed unless given. The gas viscosity defaults to that of the feed; see Gas.

source

Mocca.Framework.MoccaSystem Type

MoccaSystem <: Jutul.JutulSystem

Common supertype for every Mocca unit operation. Each unit in a flowsheet is a Jutul.SimulationModel whose system is a subtype of one of:

• DistributedUnit: has a mesh and solves PDEs (beds, membrane sides, ...)

• LumpedUnit: 0D with holdup (tanks, nodes, flash drums)

• FlowElement: no holdup, owns a flow and solves a device law (valves, flow controllers, pumps, compressors)

Implementations must provide a component_names field or override component_names.

source

Mocca.Framework.DistributedUnit Type

Unit with a mesh that solves conservation laws on cells and faces.

source

Mocca.Framework.LumpedUnit Type

0D unit with holdup.

source

Mocca.Framework.FlowElement Type

Unit without holdup that owns a flow and solves a device law.

source

Mocca.Framework.component_names Function

component_names(sys::MoccaSystem)

Names of the chemical components tracked by sys, in index order.

source

Mocca.Networks.AbstractFlowsheet Type

AbstractFlowsheet

Units and the connections between them, which setup_flowsheet_model assembles: a UnitNetwork, or a Flowsheet of the user interface, which holds one. Subtypes provide network.

source

Mocca.Networks.UnitNetwork Type

UnitNetwork()

Named units (beds, devices, ...), the connections between their ports, and the boundary conditions of device ports that are not connected. Build it with add_unit!, connect! and set_boundary!, and turn it into a Jutul.MultiModel with setup_flowsheet_model. For a process described in terms of materials and device roles, use a Flowsheet instead.

Unconnected device ports are boundaries and need a PortCondition, set with set_boundary! or per stage.

source

Mocca.Networks.network Function

network(fs) → UnitNetwork

The units and connections of fs.

source

Mocca.Networks.add_unit! Function

add_unit!(fs, name, model)

Add a unit model (a Jutul.SimulationModel of a MoccaSystem) to the flowsheet under name.

source

Mocca.Networks.connect! Function

connect!(fs, connection::AbstractConnection)
connect!(fs, (unit, port) => (device, device_port))

Add a connection to the flowsheet, after checking it with check_connection. The pair form connects port of a unit with holdup to device_port (:inlet or :outlet) of a flow device, as a StreamConnection; the pair may be given in either order.

source

Mocca.Networks.set_boundary! Function

set_boundary!(fs, device, port, condition::PortCondition)

Set the default boundary condition for an unconnected device port. Stages can override it.

source

Mocca.Networks.PortRef Type

PortRef(name, port)

Port port of the unit or device name in a flowsheet, such as PortRef(:Bed, :top). A (name, port) tuple converts to one.

source

Mocca.Networks.AbstractConnection Type

AbstractConnection

A connection between ports of the units of a flowsheet. Add one with connect!. An implementation provides:

• connection_ports(c): the PortRefs it joins;

• check_connection(net, c): throw an error if c cannot join the units of the UnitNetwork net, given the connections it already has;

• add_cross_terms!(mm, net, c): add its cross terms to the Jutul.MultiModel mm of net.

StreamConnection, the gas stream through a device port, is the one built in.

source

Mocca.Networks.connection_ports Function

The PortRefs connection c joins.

source

Mocca.Networks.check_connection Function

Throw an error if connection c cannot join the units of net.

source

Mocca.Networks.add_cross_terms! Function

Add the cross terms of connection c to the MultiModel mm built from net.

source

Mocca.Networks.StreamConnection Type

StreamConnection(unit::PortRef, device::PortRef)

The gas stream between a port of a unit with holdup and a port of a FlowDevice. Each device port takes at most one; a unit port may take several. Its cross terms are a PortStateCT on the device and a StreamCT for each of the unit's stream_equations.

source

Mocca.Networks.holdup_units Function

Units with holdup (every unit but the devices), in flowsheet order.

source

Adding a kind of unit

A new kind of unit is a MoccaSystem with these methods.

Mocca.Framework.stream_equations Function

stream_equations(model) → Tuple of Symbols

Equations of model that a stream through one of its ports adds to: the mass balance, and the energy balance if the unit solves one.

source

Mocca.Framework.setup_unit_state Function

setup_unit_state(model; pressure, temperature, composition,
    wall_temperature = temperature) → Dict

Uniform initial state of a unit with holdup. The default sets the gas state, and the wall temperature if the unit has a wall; a unit with more primary variables adds a method.

source

Mocca.Framework.setup_unit_parameters Function

setup_unit_parameters(model; overrides...) → Dict

Parameters of a unit, with overrides for any of them.

source

Mocca.Framework.holds_sorbent Function

holds_sorbent(model) → Bool

Whether model is a bed: a unit holding sorbent, which steps act on and whose sorbent counts in the productivity of a process.

source

Mocca.Balances.AxialDispersion Type

AxialDispersion

Axial dispersion coefficient [m²/s] of the gas in a distributed unit, one value for the unit. It is a parameter, so it can be changed per case, set by a design variable and differentiated without rebuilding the model.

source

Devices, schedules and steady states

Mocca.Devices.FlowDevice Type

FlowDevice(component_names)

A flow element with an inlet and an outlet port. The flow law (closed, valve, flow controller, ...) is set per stage through the law force; see AbstractDeviceLaw. Build its model with setup_device_model.

source

Mocca.Devices.setup_device_model Function

setup_device_model(component_names)

Build the SimulationModel for a FlowDevice.

source

Mocca.Devices.AbstractDeviceLaw Type

AbstractDeviceLaw

How a FlowDevice relates its flow to the states at its ports. Passed as the law force for each stage. Implementations provide device_flow.

source

Mocca.Devices.device_flow Function

device_flow(law, state, time) → F

Molar flow [mol/s] from inlet to outlet prescribed by law, given the device state (with port pressures, temperatures and compositions).

source

Mocca.Devices.LinearValve Type

LinearValve(conductance; opening = 1.0)

Valve whose volumetric flow at upstream conditions is linear in the pressure drop, q = opening·conductance·(P_in − P_out), with conductance in m³/(s·Pa). The molar flow uses the upstream concentration, so it changes sign smoothly through zero.

source

Mocca.Devices.VolumetricFlow Type

VolumetricFlow(rate)

Flow controller delivering rate [m³/s] of gas, measured at the outlet pressure and inlet temperature, from inlet to outlet. Use it for a feed specified as a superficial velocity at the bed inlet: rate = v·A.

source

Mocca.Networks.PortStateCT Type

PortStateCT(cell)

Cross term setting a device port copy to the state of unit cell. Registered with the device as target and the unit as source, on the device's :inlet_state or :outlet_state equation.

source

Mocca.Networks.StreamCT Type

StreamCT(cell, port)

Cross term adding the flow through device port port (1 = inlet, 2 = outlet) to the balances of unit cell. Registered with the unit as target and the device as source, once for the :mass_conservation equation and, if the unit has an energy balance, once for :energy_column.

The stream composition, temperature and pressure are those of the upstream device port, so component i enters or leaves at F·y_s,i and every component is conserved across the connection.

source

Mocca.Operation.Stage Type

Stage(name, duration; settings...)

One step of an operating schedule. Each keyword names a device and gives its settings for the stage as a NamedTuple with any of law, inlet and outlet (see Jutul.setup_forces for a FlowDeviceModel). Devices not named are closed. Time-dependent profiles in port conditions use time since the start of the stage.

Stage("adsorption", 15.0;
    V_feed = (law = VolumetricFlow(v_feed * A),),
    V_product = (law = LinearValve(C),))

Settings can also be passed as device => settings pairs, which is convenient when they are built programmatically: Stage("adsorption", 15.0, :V_feed => (law = Closed(),)).

source

Mocca.Operation.setup_schedule Function

setup_schedule(fs, model, stages; num_cycles = 1, max_dt = 1.0, first_dt = max_dt,
    reference_durations = nothing) → (forces, timesteps)

Forces and time steps for num_cycles repetitions of stages on the flowsheet fs and its MultiModel. Each stage is split into steps no longer than max_dt. With first_dt < max_dt, each stage starts with a step of first_dt, doubling up to max_dt, which helps the solver through the abrupt change when a valve opens or a flow controller starts. Boundary ports use the stage's port condition if given, otherwise the flowsheet default from set_boundary!.

With reference_durations, one per stage, each stage instead takes the steps of its reference duration, stretched in proportion to its own. Changing the durations then keeps the number of steps, which is how the stage duration gradients of force_gradients treat them, so an optimiser that changes durations sees the same model that the gradients describe.

source

Mocca.Operation.parallel_stages Function

parallel_stages(label => stages, ...; cycle_time = longest sequence) → Vector{Stage}

Combine stage sequences that run at the same time on separate devices into one sequence of stages. An example is two beds in series whose blowdown and evacuation times differ:

parallel_stages("silica gel" => sg_stages, "13X" => zeolite_stages)

Each sequence is padded to cycle_time with an idle stage in which its devices are closed. The combined sequence has a stage boundary wherever any sequence has one, and each of its stages holds the settings of every sequence's current stage. It is named after them, prefixed by their labels, such as "silica gel blowdown, 13X blowdown". Where a stage is split, its pressure profiles keep running from the start of the original stage.

Each device may appear in only one sequence.

source

Mocca.Devices.PortCondition Type

PortCondition(; pressure, temperature, composition)

Prescribed state at a device port that is not connected to a unit, i.e. a boundary such as a feed, a product header or a vacuum line. pressure may be a number or a time-dependent profile such as ExponentialRamp.

source

Mocca.Devices.ExponentialRamp Type

ExponentialRamp(start, stop, λ; t0 = 0.0)

Value moving from start towards stop as stop + (start − stop)·exp(−λ(t − t0)). The pressure profile used for pressurisation, blowdown and evacuation stages.

source

Mocca.Networks.port_conductance Function

port_conductance(model, cell)

Volumetric conductance [m³/(s·Pa)] of the half-cell connection between cell and the end of a distributed unit, k·A/(Δx/2)/μ. A LinearValve with this conductance reproduces the pressure boundary conditions of earlier Mocca versions.

source

Mocca.Interface.default_conductance Function

default_conductance(fs, device)

Valve conductance [m³/(s·Pa)] for device: that of the half cell at the bed end it is attached to (see port_conductance), so that the valve adds no more resistance than the bed's end and the bed follows the pressure set at the other side. For a link, the two half cells in series.

source

Mocca.Networks.setup_flowsheet_model Function

setup_flowsheet_model(fs) → Jutul.MultiModel

Assemble the flowsheet into a MultiModel, with the cross terms of every connection (see add_cross_terms!).

source

Mocca.Networks.setup_flowsheet_state Function

setup_flowsheet_state(fs; kwarg...) → Dict

Initial state for every unit. Pass the initial state of each unit with holdup as a keyword named after the unit, e.g. Bed = setup_unit_state(fs[:Bed]; ...). Device states are initialised from what their ports are attached to.

source

Mocca.Networks.setup_flowsheet_parameters Function

setup_flowsheet_parameters(fs; kwarg...) → Dict

Parameters for every unit, from setup_unit_parameters; pass a Dict or NamedTuple of overrides per unit as a keyword named after it.

source

Mocca.Physics.DualSiteLangmuir Type

DualSiteLangmuir{N, T} <: AbstractIsotherm

Dual-site Langmuir isotherm (Haghpanah et al. 2013), concentration basis.

q*ᵢ = qsbᵢ·bᵢ·Cᵢ/(1 + Σⱼ bⱼ·Cⱼ) + qsdᵢ·dᵢ·Cᵢ/(1 + Σⱼ dⱼ·Cⱼ)

Two adsorption sites per component:

• Site b (high energy): saturation qsb, affinity b(T) = b0·exp(-ΔUb/(R·T))

• Site d (low energy): saturation qsd, affinity d(T) = d0·exp(-ΔUd/(R·T))

Enthalpy is state-independent, precomputed for each component as the average of its two site energies weighted by their saturation capacities.

source

Mocca.Physics.AbstractIsotherm Type

AbstractIsotherm

Abstract base type for adsorption isotherm models.

Implementations must provide:

• compute_equilibrium(iso, C, T): equilibrium loading per component [mol/m³ adsorbent]

• compute_enthalpy(iso, C, T): isosteric heat of adsorption per component [J/mol]

source

Mocca.Physics.AbstractMassTransfer Type

AbstractMassTransfer

Abstract base type for mass transfer models.

Implementations must provide:

• compute_mass_transfer_rate(mt, C, q, qstar): Calculate mass transfer rate

source

Mocca.Physics.LinearDrivingForce Type

LinearDrivingForce{T, D} <: AbstractMassTransfer

Linear Driving Force (LDF) mass transfer model with single lumped resistance.

The rate equation is:

∂q/∂t = k * (q* - q)

where k = 15 * ε_p * D_p / r_p² * C / q* and D_p = D_m / τ. A component the sorbent does not take up (q* = 0 at any concentration, such as one without isotherm data) has no mass transfer.

source

Mocca.Operation.simulate_to_cyclic_steady_state Function

simulate_to_cyclic_steady_state(model, state0, parameters, forces, timesteps;
    tol = 1e-4, max_cycles = 100, acceleration = :anderson, depth = 5,
    progress = false, info_level = -1, kwarg...)

Repeat one cycle, given by its forces and timesteps, until the cycle_change of a cycle falls below tol or max_cycles is reached. Works for a single unit model and for flowsheets.

Each cycle maps its initial state x to its final state Φ(x), and cyclic steady state is the fixed point x = Φ(x). With acceleration = :none the next cycle starts from Φ(x) (successive substitution). With :anderson it starts from a combination of the last depth cycles that extrapolates towards the fixed point (Anderson acceleration), which usually needs far fewer cycles. The extrapolated state is kept within the bounds of each variable, and if a cycle started from it fails, the next cycle starts from Φ(x) instead.

With progress (true, an IO or a CycleProgress), each cycle is printed with its change and how far it has come towards tol.

Returns a NamedTuple with the states and time steps of the last cycle, its initial state state0, the change per cycle in history, the number of cycles simulated, and whether it converged. Extra keywords go to simulate_process.

source

Mocca.Operation.newton_cyclic_steady_state Function

newton_cyclic_steady_state(model, state0, parameters, forces, timesteps;
    tol = 1e-10, max_iterations = 10, krylov = :auto, growth = 100, max_backtracks = 6,
    keep_steps_below = 1e-3, progress = false, info_level = -1, kwarg...)

Converge a cyclic steady state with Newton's method from state0, such as the end of a few plain cycles. The cycle, given by its forces and timesteps, maps its initial state x to its final state Φ(x), and each iteration solves

(I − ∂Φ/∂x) δ = Φ(x) − x,   x ← x + δ.

Products with ∂Φ/∂x are forward sensitivities over the cycle as simulated: a change in the start carried through the linearised steps. Each step is assembled and factorised once, so that many products cost little more than one. Convergence is quadratic, so a few iterations reach a steady state far tighter than simulate_to_cyclic_steady_state can, which needs many cycles when the slowest mode of the cycle decays slowly. A rough start is enough, once a few plain cycles have let the fastest modes die out: the control of the steps described below does the rest. Start each stage with short steps (first_dt in setup_schedule) so that the cycle map is smooth; see cyclic_steady_state_gradient.

With krylov = false the Jacobian ∂Φ/∂x is built in full, one column for each primary variable of the units with holdup, in one sweep over the cycle. With jacobian, the jacobian returned by an earlier call (for example for slightly different settings, as in an optimisation), the iterations start from it instead of building one, and improve it with a Broyden update from each iteration's step and change in residual (a quasi-Newton method), each costing one cycle. The Jacobian is rebuilt whenever an iteration reduces the cycle change by less than a factor 1/refresh_ratio.

With krylov = true the Jacobian is never built: each iteration solves for the step with GMRES (Jacobian-free Newton–Krylov), with one forward sensitivity per GMRES iteration, to a relative residual that tightens with the cycle change. GMRES needs about as many products as the cycle has modes that do not die out within a cycle, so this is cheaper than the full Jacobian when there are many more variables than slow modes. The products also give an approximation of the Jacobian on the directions GMRES explored, which preconditions the later iterations and is returned as jacobian; a jacobian given preconditions the first. At most max_products products are taken per iteration (by default, the number of variables). With krylov = :auto (the default), GMRES is used when the units with holdup have more than 800 primary variables (a bed of about 130 cells with two components, 100 with three), and the full Jacobian otherwise; krylov = n sets that number of variables to n instead.

Far from the steady state the full step can overshoot: the slow modes of the cycle make it long. Each step is tried, and taken if the cycle simulates from it and the norm of the residual Φ(x) − x (in the units of the cycle change) is less than growth times that of the best iterate so far. A Newton step often raises the residual of a few cells before the next one converges, so only steps that make it much worse are cut back: a quasi-Newton step first gets a new Jacobian, then the step is shortened, from a quadratic fit of the residual along it, or to a quarter if the cycle failed to simulate, up to max_backtracks times. If no step is accepted the next iteration starts from the end of the cycle, as a plain cycle would. The cycle simulated for an accepted step is that of the next iteration, so a full step costs no extra cycle.

Where the solver has to cut a time step can depend on where the cycle starts, and a cut that comes and goes between nearby starts makes the cycle map jump by about the time-stepping error, which near the steady state is far more than the change left (on the wet LPP VSA, a few steps at the start of pressurisation from vacuum). Once the cycle change is below keep_steps_below, every later cycle therefore takes the steps of that iteration's cycle, cuts included (and adds any further cut it needs), so that the map is smooth. The steady state found is that of those steps, which differ from the steps given only where the solver cut them.

Stops when the cycle_change is below tol. With progress (true, an IO or a CycleProgress), each iteration is printed with its cycle change and its step if shortened, and each Jacobian built, GMRES solve or step rejected. Returns the same fields as simulate_to_cyclic_steady_state, with cycles counting the iterations, plus jacobian, the last Jacobian used (as the inverse of I − ∂Φ/∂x), jacobian_builds, how many were built, products, the number of GMRES products, trials, the number of steps rejected, and fallbacks, the iterations that fell back to a plain cycle. Extra keywords go to setup_process_simulator.

source

Mocca.Operation.simulate_until_mass_balance Function

simulate_until_mass_balance(model, state0, parameters, forces, timesteps;
    inflow, outflow, component = 1, tol = 0.005, consecutive = 5,
    min_cycles = 50, max_cycles = 300, progress = false, info_level = -1, kwarg...)

Repeat one cycle of a flowsheet, given by its forces and timesteps, until the mass balance of component over the flowsheet closes: the relative difference between the moles entering through the inflow devices and leaving through the outflow devices is below tol for consecutive cycles in a row, after at least min_cycles. This is the cyclic steady state criterion of Haghpanah et al. (2013) and Krishnamurthy et al. (2014), with CO2 as the component.

It only checks one component. A more strongly adsorbed component, such as water, can still be accumulating when it is met; use simulate_to_cyclic_steady_state to require every variable to settle.

With progress (true, an IO or a CycleProgress), each cycle is printed with its balance error, how many cycles in a row have been below tol, and an estimate of the cycles left.

Returns a NamedTuple with the states and time steps of the last cycle, its initial state state0, the balance error of each cycle in history, the number of cycles simulated, and whether the criterion was met (converged). Extra keywords go to simulate_process.

source

Mocca.Operation.CycleProgress Type

CycleProgress(io = stdout; live = io isa Base.TTY)

Printout of the progress of the cyclic steady state solvers, given to them as progress. Each solver is a phase with a tolerance of its own (such as Anderson acceleration to 1e-3, then Newton's method to the final tolerance), and each phase has a line with a bar. The line shows the cycle count, the measure of change against the phase's tolerance, a bar of how far the phase has come (on a log scale, from its first cycle to its tolerance, by the smallest change so far, as the change of an accelerated cycle goes up and down), an estimate of the cycles left from the recent rate of decrease of that smallest change, and the time taken. When a phase ends, its line says whether it converged, or for a phase that hands over to the next, that it is done.

With live, as in a terminal, the lines of the phases planned with plan_phases! are drawn together and rewritten in place, the later ones waiting until their turn. Otherwise each cycle gets a line of its own.

source

Mocca.Operation.cycle_change Function

cycle_change(model, state_start, state_end) → Float64

Largest relative change of any primary variable of a unit with holdup between the start and the end of a cycle: max |x_end − x_start| / max |x_start|, taken per variable, and per component for mole fractions and loadings, so that a trace component is measured against its own level. A mole fraction is measured against at least 10⁻⁶, and a loading against at least 0.1% of the component's saturation capacity, so that an absent component does not measure round-off against itself. Devices are skipped, since they hold no inventory.

source

Mocca.Operation.cyclic_steady_state_gradient Function

cyclic_steady_state_gradient(model, state0, parameters, forces, timesteps, G;
    targets = nothing, with_forces = false, rtol = 1e-8, maxiter = 100,
    nonlinear_tolerance = 1e-8, kwarg...) → NamedTuple

Gradient at cyclic steady state of the objective G over one cycle, given by its forces and timesteps. state0 is the state at the start of a cycle at steady state, such as state0 from newton_cyclic_steady_state or simulate_to_cyclic_steady_state. G is a Jutul sum objective, such as a StreamObjective, or a NamedTuple of them, which share the simulation of the cycle.

The cycle maps its initial state x to its final state Φ(x, p), and the steady state x* solves x* = Φ(x*, p), so it moves with any parameter or setting p. By the implicit function theorem, the gradient of J(p) = G(x*(p), p) is

dJ/dp = ∂J/∂p + μᵀ ∂Φ/∂p,   where   (I − ∂Φ/∂x)ᵀ μ = ∂J/∂x.

Each product with (∂Φ/∂x)ᵀ is an adjoint solve over the cycle, and GMRES finds μ from them; rtol and maxiter control it. Pass the jacobian from newton_cyclic_steady_state to precondition it, which cuts the iterations from tens to a few. The gradient is only as accurate as the steady state: an error e in the cycle change gives a gradient error of about e/(1 − ρ), where ρ is the decay factor per cycle of the slowest mode. newton_cyclic_steady_state converges it far enough.

The cycle map must be smooth for the gradient to mean anything. Abrupt stage starts taken in one long step can make it jump: the solver then lands on a different solution of that step for a slightly different start. Start each stage with short steps (first_dt in setup_schedule).

Returns a NamedTuple (one per objective, if G is a NamedTuple) with:

• objective: the value of G over the cycle;

• parameters: the gradient with respect to the parameters of each unit, laid out like the output of Jutul.solve_adjoint_sensitivities, for the parameters in targets (a Dict of unit => parameter names) or all of them;

• forces and durations: if with_forces, the gradient with respect to each stage's device settings and duration, as in force_gradients. A longer stage also lengthens the cycle; divide by the cycle time for rates;

• state0: the multiplier μ, laid out like parameters: the change in the objective at steady state per unit of each initial-state variable added at the start of a cycle;

• iterations and converged for GMRES.

Extra keywords go to setup_process_simulator.

source

Mocca.Operation.force_gradients Function

force_gradients(model, state0, parameters, forces, timesteps, G;
    nonlinear_tolerance = 1e-8, kwarg...) → NamedTuple

Gradient of the objective G with respect to the device settings and the duration of each stage, from one simulation from state0. G is a Jutul sum objective, such as a StreamObjective.

Returns the objective value and, with one entry per stage (per distinct force set, in order of use):

• forces: each entry has the layout of the stage's forces, with every number replaced by the derivative of the objective with respect to it: g.forces[2][:V_feed].law.rate is the derivative with respect to the feed rate in the second stage. Mole fractions are treated as independent.

• durations: the derivative of the objective with respect to the duration of the stage [per s]. The stage is stretched with its steps in proportion; pressure ramps run on the time since the start of their stage, so the later stages are the same in their own time. This is the derivative of the simulation as discretised, the change seen with reference_durations in setup_schedule; without it, a new duration keeps the stage's first steps, which differs by the time-stepping error. An objective that depends on the time other than through the time step is differentiated as if it did not.

Extra keywords go to setup_process_simulator.

source

Mocca.Operation.StreamObjective Type

StreamObjective(device; component = 1, scale = 1.0)

Objective for adjoint gradients: scale times the moles of component that cross device from its inlet to its outlet, summed over the steps as in stream_totals.

source

Mocca.Operation.VacuumPumpObjective Type

VacuumPumpObjective(device; discharge_pressure = 101325.0, efficiency = 0.72,
    γ = 1.4, scale = 1.0)

Objective for adjoint gradients: scale times the vacuum pump work [J] for the gas leaving through device, as in vacuum_pump_energy.

source

Mocca.Operation.CompressorObjective Type

CompressorObjective(device; suction_pressure = 101325.0, efficiency = 0.72, γ = 1.4, scale = 1.0)

Objective for adjoint gradients: scale times the work [J] to compress the gas entering through device from suction_pressure up to the pressure at its outlet, such as a feed compressor taking flue gas at atmospheric pressure up to the bed inlet pressure, as in compressor_energy. Only steps with forward flow and an outlet above the suction pressure count. The gas is compressed adiabatically with isentropic efficiency, from the device's inlet temperature.

source

Mocca.Framework.SimulationFailure Type

SimulationFailure(message)

The solver could not complete a simulation: Jutul gave up on a time step after cutting it, so the states stop short of the end of the schedule.

source