3.4. wavesim.iteration

class Domain(pixel_size, shape, device)[source]

Bases: object

Base class for all simulation domains

This base class defines the interface for operations that are common for all simulation types, and for MultiDomain. todo: the design using slots minimizes memory use, but it is a suboptimal design because it mixes mutable and immutable state. This design should be revisited so that the Domain is immutable, and the code that runs the algorithms performs the memory management.

abstractmethod add_source(slot, weight)[source]
abstractmethod clear(slot)[source]

Clears the data in the specified slot

coordinates(dim, type='linear')[source]

Returns the real-space coordinates along the specified dimension, starting at 0

coordinates_f(dim)[source]

Returns the Fourier-space coordinates along the specified dimension

abstractmethod create_empty_vdot()[source]

Create an empty tensor for the Vdot tensor

abstractmethod get(slot, copy=False)[source]

Returns the data in the specified slot.

Parameters:
  • slot (int) – slot from which to return the data

  • copy – if True, returns a copy of the data. Otherwise, may return the original data possible.

Note that this data may be overwritten by the next call to domain.

abstractmethod inner_product(slot_a, slot_b)[source]

Computes the inner product of two data vectors

Note

The vectors may be represented as multidimensional arrays, but these arrays must be contiguous for this operation to work. Although it would be possible to use flatten(), this would create a copy when the array is not contiguous, causing a hidden performance hit.

abstractmethod inverse_propagator(slot_in, slot_out)[source]

Applies the operator (L+1) x .

This operation is not needed for the Wavesim algorithm, but is provided for testing purposes, and can be used to evaluate the residue of the solution.

abstractmethod medium(slot_in, slot_out, mnum)[source]

Applies the operator 1-Vscat.

abstractmethod mix(weight_a, slot_a, weight_b, slot_b, slot_out)[source]

Mixes two data arrays and stores the result in the specified slot

abstractmethod propagator(slot_in, slot_out)[source]

Applies the operator (L+1)^-1 x.

abstractmethod set(slot, data)[source]

Copy the date into the specified slot

abstractmethod set_source(source)[source]

Sets the source term for this domain.

domain_operator(domain, function, **kwargs)[source]

Constructs various operators by combining calls to ‘medium’, ‘propagator’, etc.

todo: this is currently very inefficient because of the overhead of copying data to and from the device

forward(domain, slot_in, slot_out)[source]

Evaluates the forward operator A= c⁻¹ (L + V)

Parameters:
  • domain (Domain) – Domain object

  • slot_in (int) – slot holding input x. This slot will be overwritten!

  • slot_out (int) – output slot that will receive A x

is_zero(x)[source]

Check if x is zero

Some functions allow specifying 0 or 0.0 instead of a torch tensor, to indicate that the array should be cleared. This function returns True if x is a scalar 0 or 0.0. It raises an error if x is a scalar that is not equal to 0 or 0.0, and returns False otherwise.

preconditioned_iteration(domain, slot_in=0, slot_out=0, slot_tmp=1, alpha=0.75, compute_norm2=False)[source]

Run one preconditioned iteration.

Parameters:
  • domain – Domain object

  • slot_in (int) – slot holding input x. This slot will be overwritten!

  • slot_out (int) – output slot that will receive the result

  • slot_tmp (int) – slot for temporary storage. Cannot be equal to slot_in, may be equal to slot_out

  • alpha – relaxation parameter for the Richardson iteration

  • compute_norm2 – when True, returns the squared norm of the residual. Otherwise, returns 0.0

Richardson iteration:

x -> x + α (y - A x)

Preconditioned Richardson iteration:
x -> x + α Γ⁻¹ (y - A x)

= x + α c B (L+1)⁻¹ (y - A x) = x + α c B (L+1)⁻¹ (y - c⁻¹ [L+V] x) = x + α c B (L+1)⁻¹ (y + c⁻¹ [1-V] x - c⁻¹ [L+1] x) = x + α B [(L+1)⁻¹ (c y + B x) - x] = x - α B x + α B (L+1)⁻¹ (c y + B x)

preconditioned_operator(domain, slot_in, slot_out)[source]

Evaluates the preconditioned operator Γ⁻¹ A

Where Γ⁻¹ = c B (L+1)⁻¹

Note: the scale factor c that makes A accretive and V a contraction is

included in the preconditioner. The Richardson step size is _not_.

Operator A is the original non-scaled operator, and we have (L+V) = c A Then:

Γ⁻¹ A = c B (L+1)⁻¹ A

= c B (L+1)⁻¹ c⁻¹ (L+V) = B (L+1)⁻¹ (L+V) = B (L+1)⁻¹ ([L+1] - [1-V]) = B - B (L+1)⁻¹ B

Parameters:
  • domain (Domain) – Domain object

  • slot_in (int) – slot holding input x. This slot will be overwritten!

  • slot_out (int) – output slot that will receive A x

preconditioner(domain, slot_in, slot_out)[source]

Evaluates Γ⁻¹ = c B(L+1)⁻¹

Parameters:
  • domain (Domain) – Domain object

  • slot_in (int) – slot holding input x. This slot will be overwritten!

  • slot_out (int) – output slot that will receive A x

run_algorithm(domain, source, alpha=0.75, max_iterations=1000, threshold=1.e-6, full_residuals=False)[source]

WaveSim update

Parameters:
  • domain (Domain) – Helmholtz base parameters

  • source – source field

  • alpha – relaxation parameter for the Richardson iteration

  • max_iterations – maximum number of iterations

  • threshold – threshold for the residual norm

  • full_residuals – when True, returns list of residuals for all iterations. Otherwise, returns final residual

Returns:

u, iteration count, residuals