darsia.BeckmannProblem#

class BeckmannProblem(grid, weight=None, options={})[source]#

Bases: EMD

Base class for setting up the Beckmann problem.

The Beckmann problem is the optimisation problem

\[\inf \|u\|_{L^1} \quad \text{s.t.} \quad \operatorname{div} u = m_2 - m_1, \; u \in H(\operatorname{div}),\]

where \(u\) is the flux and \(m_1\), \(m_2\) are the mass distributions transported from \(m_1\) to \(m_2\). It is solved approximately with an iterative TPFA-type finite-volume method, exploiting the connection to the lowest-order Raviart-Thomas mixed finite element.

Two solution strategies are implemented in specialised subclasses:

a finite-volume quasi-Newton method (BeckmannNewtonSolver);
a finite-volume split-Bregman method (BeckmannBregmanSolver).

Initialise the variational Wasserstein distance.

Parameters:

griddarsia.Grid

Tensor grid associated with the images.

weightdarsia.Image, optional

Per-pixel weight.

optionsdict, optional

Options for the nonlinear and linear solver. Recognised keys, with their defaults:

l1_mode (darsia.L1Mode)      "RAVIART_THOMAS"
    "RAVIART_THOMAS"              exact RT0 integration (MFEM)
    "CONSTANT_SUBCELL_PROJECTION" subcell L2 projection (mixed FV)
    "CONSTANT_CELL_PROJECTION"    cell-based L2 projection
mobility_mode (MobilityMode) MobilityMode.CELL_BASED
    CELL_BASED / CELL_BASED_ARITHMETIC / CELL_BASED_HARMONIC /
    SUBCELL_BASED / FACE_BASED    face flux-norm averaging scheme
num_iter (int)               100     max iterations
tol_residual (float)         inf     residual tolerance
tol_increment (float)        inf     increment tolerance
tol_distance (float)         inf     distance tolerance
L (float)                    1.0     regularisation / approx flux norm
linear_solver (str)          "direct"  ("direct"|"amg"|"cg"|"ksp")
formulation (str)            "pressure"
    ("full"|"flux_reduced"|"pressure")   what is eliminated
linear_solver_options (dict) {}
amg_options (dict)           {}
aa_depth (int)               0       Anderson acceleration depth
aa_restart (int)             None    Anderson acceleration restart
regularization (float)       eps     avoid division by zero
lumping (bool)               True    lump the mass matrix

Methods

`broken_darcy_with_custom_flux_block`(...)	Construct the broken Darcy operator with given flux-flux block.
`cell_weighted_flux`(cell_flux)	Compute the cell-weighted flux.
`distance_matrix`(images)	Compute the distance between each iteam of a list.
`eliminate_flux`(jacobian, residual)	Eliminate the flux block from the jacobian and residual.
`eliminate_lagrange_multiplier`(...)	Eliminate the lagrange multiplier from the reduced system.
`exact_linearization`(solution)	Compute the exact linearization of the constrained minimization problem.
`flat_view`(img)	Flatten the image to a vector.
`flux_view`(vector)	Extract the flux from the vector.
`l1_dissipation`(flat_flux)	Compute the l1 dissipation of the solution.
`linear_solve`(matrix, rhs[, ...])	Solve the linear system.
`ndofs`()	Return the total number of degrees of freedom.
`optimality_conditions`(solution, ...)	Evaluate optimality conditions of the constrained minimization problem.
`pressure_view`(vector)	Extract the pressure from the vector.
`rescaled_flux_optimality_conditions`(solution)	Evaluate scaled optimality conditions of the constrained minimization problem.
`solve_beckmann_problem`(mass_diff)	Solve for the Wasserstein distance.
`transport_density`(flat_flux[, weighted, flatten])	Compute the transport density from the solution.
`transport_density_faces`(flat_flux)	Compute the transport density from the solution.

Attributes

`grid`	darsia.Grid: grid
`voxel_size`	np.ndarray: voxel size
`options`	dict: options for the solver
`regularization`	float: regularization parameter
`verbose`	bool: verbosity
`mobility_mode`	str: mode for computing the mobility
`weight`	Weight defined on cells
`callbacks`	list: list of callbacks to be called during the optimization

broken_darcy_with_custom_flux_block(flux_flux_block)[source]#

Construct the broken Darcy operator with given flux-flux block.

Parameters:

flux_flux_blocksps.csc_matrix: Flux-flux block.

Returns:

sps.csc_matrix: Broken Darcy operator with given flux-flux block.

cell_weighted_flux(cell_flux)[source]#

Compute the cell-weighted flux.

Parameters:

cell_fluxnp.ndarray: Cell fluxes.

Returns:

np.ndarray: Cell-weighted flux.

eliminate_flux(jacobian, residual)[source]#

Eliminate the flux block from the jacobian and residual.

Employ a Schur complement/block Gauss elimination approach.

Parameters:

jacobiansps.csc_matrix: Jacobian.
residualnp.ndarray: Residual.

Returns:

tuple: Reduced jacobian, reduced residual, inverse of flux block.

eliminate_lagrange_multiplier(reduced_jacobian, reduced_residual)[source]#

Eliminate the lagrange multiplier from the reduced system.

Employ a Schur complement/block Gauss elimination approach.

Parameters:

reduced_jacobiansps.csc_matrix: Reduced jacobian.
reduced_residualnp.ndarray: Reduced residual.

Returns:

tuple: Fully reduced jacobian, fully reduced residual.

exact_linearization(solution)[source]#

Compute the exact linearization of the constrained minimization problem.

Parameters:

solutionnp.ndarray: Solution.

Returns:

sps.csc_matrix: Exact linearization.

flat_view(img)[source]#

Flatten the image to a vector.

Parameters:

imgnp.ndarray: Image.

Returns:

np.ndarray: Flattened image.

flux_view(vector)[source]#

Extract the flux from the vector.

Parameters:

vectornp.ndarray: Vector.

Returns:

np.ndarray: Flux.

l1_dissipation(flat_flux)[source]#

Compute the l1 dissipation of the solution.

Parameters:

flat_fluxnp.ndarray: Flat fluxes.

Returns:

float: L1 dissipation.

linear_solve(matrix, rhs, previous_solution=None, reuse_solver=False)[source]#

Solve the linear system.

Defines the Schur complement reduction and the pure pressure reduction, if selected. For reusing the setup, the resulting solver is cached as self.linear_solver.

Parameters:

matrixsps.csc_matrix: Matrix.
rhsnp.ndarray: Right hand side.
previous_solutionnp.ndarray: Previous solution. Defaults to None.
reuse_solverbool: Reuse the solver. Defaults to False.

Returns:

tuple: Solution, stats.

ndofs()[source]#

Return the total number of degrees of freedom.

Returns:

int: Total number of degrees of freedom.

optimality_conditions(solution, beckmann_problem_rhs)[source]#

Evaluate optimality conditions of the constrained minimization problem.

Parameters:

solutionnp.ndarray: Solution.
beckmann_problem_rhsnp.ndarray: Right hand side of the Beckmann problem.

Returns:

np.ndarray: Residual.

pressure_view(vector)[source]#

Extract the pressure from the vector.

Parameters:

vectornp.ndarray: Vector.

Returns:

np.ndarray: Pressure.

rescaled_flux_optimality_conditions(solution)[source]#

Evaluate scaled optimality conditions of the constrained minimization problem.

Scale the flux equation by the face weights divided by the total distance. This results in no division by zero.

abstractmethod solve_beckmann_problem(mass_diff)[source]#

Solve for the Wasserstein distance.

Parameters:

mass_diffnp.ndarray: Difference of the two distributions.

Returns:

tuple: Distance, solution, info.

transport_density(flat_flux, weighted=True, flatten=True)[source]#

Compute the transport density from the solution.

Parameters:

flat_fluxnp.ndarray: Face fluxes.
weightedbool: Apply weighting. Defaults to True.
flattenbool: Flatten the result. Defaults to True.

Returns:

np.ndarray: Transport density, flattened if requested.

transport_density_faces(flat_flux)[source]#

Compute the transport density from the solution.

Parameters:

flat_fluxnp.ndarray: Face fluxes.

Returns:

np.ndarray: Transport density.

callbacks#: list: list of callbacks to be called during the optimization

grid#: darsia.Grid: grid

mobility_mode#: str: mode for computing the mobility

options#: dict: options for the solver

regularization#: float: regularization parameter

verbose#: bool: verbosity

voxel_size#: np.ndarray: voxel size

weight#: Weight defined on cells