Operators

Operators#

MOLE provides a comprehensive set of mimetic operators for numerical computations. These operators are designed to preserve important mathematical properties of the continuous operators they approximate.

Gradient Operator#

The Gradient operator computes the gradient of a scalar field in the MOLE library.

API Reference#

class Gradient : public sp_mat#

Mimetic Gradient operator.

Supports both non-periodic and periodic boundary conditions. The BC-aware constructors accept Robin coefficient vectors dc and nc representing a0 and b0 in the condition a0*U + b0*dU/dn = g. An axis is treated as periodic when all of its dc and nc entries are zero.

Public Functions

Gradient(u16 k, u32 m, Real dx)#

1-D Mimetic Gradient (non-periodic)

@argument k Order of accuracy @argument m Number of cells @argument dx Spacing between cells

Gradient(u16 k, u32 m, u32 n, Real dx, Real dy)#

2-D Mimetic Gradient (non-periodic)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction

Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz)#

3-D Mimetic Gradient (non-periodic)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument o Number of cells in z-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction @argument dz Spacing between cells in z-direction

Gradient(u16 k, u32 m, Real dx, const ivec &dc, const ivec &nc)#

1-D Mimetic Gradient (periodic or non-periodic)

@argument k Order of accuracy @argument m Number of cells @argument dx Spacing between cells @argument dc Robin coefficient a0; 2-element vector [left, right]. All-zero → periodic BC. @argument nc Robin coefficient b0; 2-element vector [left, right]. All-zero → periodic BC.

Periodic result is m×m; non-periodic result is (m+1)×(m+2).

Gradient(u16 k, u32 m, u32 n, Real dx, Real dy, const ivec &dc, const ivec &nc)#

2-D Mimetic Gradient (periodic or non-periodic per axis)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction @argument dc Robin coefficient a0; 4-element integer vector [left, right, bottom, top]. Entries 0-1 all-zero → periodic in x. Entries 2-3 all-zero → periodic in y. @argument nc Robin coefficient b0; 4-element integer vector [left, right, bottom, top].

Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz, const ivec &dc, const ivec &nc)#

3-D Mimetic Gradient (periodic or non-periodic per axis)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument o Number of cells in z-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction @argument dz Spacing between cells in z-direction @argument dc Robin coefficient a0; 6-element integer vector [left, right, bottom, top, front, back]. Entries 0-1 all-zero → periodic in x. Entries 2-3 all-zero → periodic in y. Entries 4-5 all-zero → periodic in z. @argument nc Robin coefficient b0; 6-element integer vector, same ordering as dc.

vec getP()#

Returns the weights used in the Mimetic Gradient Operators.

Note

For informational purposes only; meaningful for non-periodic 1-D gradients. Returns an empty vector for periodic or multi-D cases.

Divergence Operator#

The Divergence operator computes the divergence of a vector field in the MOLE library.

API Reference#

class Divergence : public sp_mat#

Mimetic Divergence operator.

Supports both non-periodic and periodic boundary conditions. The BC-aware constructors accept Robin coefficient vectors dc and nc representing a0 and b0 in the condition a0*U + b0*dU/dn = g. An axis is treated as periodic when all of its dc and nc entries are zero. WARNING: At the 8th order, the weight matrix Q loses positive definiteness, so the inner product induced by Q is no longer well-defined. If the inner product is not valid, the discrete integration by parts identity has no meaning, which breaks the structure that makes the divergence mimetic.

Public Functions

Divergence(u16 k, u32 m, Real dx)#

1-D Mimetic Divergence (non-periodic)

@argument k Order of accuracy @argument m Number of cells @argument dx Spacing between cells

Divergence(u16 k, u32 m, u32 n, Real dx, Real dy)#

2-D Mimetic Divergence (non-periodic)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction

Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz)#

3-D Mimetic Divergence (non-periodic)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument o Number of cells in z-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction @argument dz Spacing between cells in z-direction

Divergence(u16 k, u32 m, Real dx, const ivec &dc, const ivec &nc)#

1-D Mimetic Divergence (periodic or non-periodic)

@argument k Order of accuracy @argument m Number of cells @argument dx Spacing between cells @argument dc Robin coefficient a0; 2-element integer vector [left, right]. All-zero → periodic BC. @argument nc Robin coefficient b0; 2-element integer vector [left, right]. All-zero → periodic BC.

Periodic result is m×m; non-periodic result is (m+2)×(m+1).

Divergence(u16 k, u32 m, u32 n, Real dx, Real dy, const ivec &dc, const ivec &nc)#

2-D Mimetic Divergence (periodic or non-periodic per axis)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction @argument dc Robin coefficient a0; 4-element integer vector [left, right, bottom, top]. Entries 0-1 all-zero → periodic in x. Entries 2-3 all-zero → periodic in y. @argument nc Robin coefficient b0; 4-element integer vector [left, right, bottom, top].

Divergence(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz, const ivec &dc, const ivec &nc)#

3-D Mimetic Divergence (periodic or non-periodic per axis)

@argument k Order of accuracy @argument m Number of cells in x-direction @argument n Number of cells in y-direction @argument o Number of cells in z-direction @argument dx Spacing between cells in x-direction @argument dy Spacing between cells in y-direction @argument dz Spacing between cells in z-direction @argument dc Robin coefficient a0; 6-element integer vector [left, right, bottom, top, front, back]. Entries 0-1 all-zero → periodic in x. Entries 2-3 all-zero → periodic in y. Entries 4-5 all-zero → periodic in z. @argument nc Robin coefficient b0; 6-element integer vector, same ordering as dc.

vec getQ()#

Returns the weights used in the Mimetic Divergence Operators.

Note

For informational purposes only; meaningful for non-periodic 1-D divergence operators. Returns an empty vector for periodic or multi-D cases.

Laplacian Operator#

The Laplacian operator computes the Laplacian of a scalar field in the MOLE library.

API Reference#

class Laplacian : public sp_mat#

Mimetic Laplacian operator.

Public Functions

Laplacian(u16 k, u32 m, Real dx)#

1-D Mimetic Laplacian Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells

  • dx – Spacing between cells

Laplacian(u16 k, u32 m, u32 n, Real dx, Real dy)#

2-D Mimetic Laplacian Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells in x-direction

  • n – Number of cells in y-direction

  • dx – Spacing between cells in x-direction

  • dy – Spacing between cells in y-direction

Laplacian(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz)#

3-D Mimetic Laplacian Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells in x-direction

  • n – Number of cells in y-direction

  • o – Number of cells in z-direction

  • dx – Spacing between cells in x-direction

  • dy – Spacing between cells in y-direction

  • dz – Spacing between cells in z-direction

Interpolation Operator#

The Interpol class performs interpolation operations in the MOLE library.

API Reference#

class Interpol : public sp_mat#

Mimetic Interpolator operator.

Public Functions

Interpol(u32 m, Real c)#

1-D Mimetic Interpolator Constructor

Parameters:

  • m – Number of cells

  • c – Weight for ends, can be any value from 0.0<=c<=1.0

Interpol(u32 m, u32 n, Real c1, Real c2)#

2-D Mimetic Interpolator Constructor

Parameters:

  • m – Number of cells in x-direction

  • n – Number of cells in y-direction

  • c1 – Weight for ends in x-direction, can be any value from 0.0<=c<=1.0

  • c2 – Weight for ends in y-direction, can be any value from 0.0<=c<=1.0

Interpol(u32 m, u32 n, u32 o, Real c1, Real c2, Real c3)#

3-D Mimetic Interpolator Constructor

Parameters:

  • m – Number of cells in x-direction

  • n – Number of cells in y-direction

  • o – Number of cells in z-direction

  • c1 – Weight for ends in x-direction, can be any value from 0.0<=c<=1.0

  • c2 – Weight for ends in y-direction, can be any value from 0.0<=c<=1.0

  • c3 – Weight for ends in z-direction, can be any value from 0.0<=c<=1.0

Interpol(bool type, u32 m, Real c)#

1-D Mimetic Interpolator Constructor

Parameters:

  • type – Dummy holder to trigger overloaded function

  • m – Number of cells

  • c – Weight for ends, can be any value from 0.0<=c<=1.0

Interpol(bool type, u32 m, u32 n, Real c1, Real c2)#

2-D Mimetic Interpolator Constructor

Parameters:

  • type – Dummy holder to trigger overloaded function

  • m – Number of cells in x-direction

  • n – Number of cells in y-direction

  • c1 – Weight for ends in x-direction, can be any value from 0.0<=c<=1.0

  • c2 – Weight for ends in y-direction, can be any value from 0.0<=c<=1.0

Interpol(bool type, u32 m, u32 n, u32 o, Real c1, Real c2, Real c3)#

3-D Mimetic Interpolator Constructor

Parameters:

  • type – Dummy holder to trigger overloaded function

  • m – Number of cells in x-direction

  • n – Number of cells in y-direction

  • o – Number of cells in z-direction

  • c1 – Weight for ends in x-direction, can be any value from 0.0<=c<=1.0

  • c2 – Weight for ends in y-direction, can be any value from 0.0<=c<=1.0

  • c3 – Weight for ends in z-direction, can be any value from 0.0<=c<=1.0

Usage Examples#

Transport Example (Gradient & Divergence)#

Transport 1D Example using Gradient and Divergence (examples/cpp/transport1D.cpp)#
 1/**
 2 * This example uses MOLE to solve the 1D advection-reaction-dispersion
 3 * equation:
 4 * https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-22.html
 5 */
 6
 7#include "mole.h"
 8#include <iostream>
 9
10int main() {
11
12  int k = 2;             // Operators' order of accuracy
13  Real a = 0;            // Left boundary
14  Real b = 130;          // Right boundary
15  int m = 26;            // Number of cells
16  Real dx = (b - a) / m; // Cell's width [m]
17  Real t = 4;            // Simulation time [years]
18  int iter = 208;        // Number of iterations
19  Real dt = t / iter;    // Time step
20  Real dis = 5;          // Dispersivity [m]
21  Real vel = 15;         // Pore-water flow velocity [m/year]
22  Real R = 2.5;          // Retardation (Cl^- = 1, Na^+ = 2.5)
23  Real C0 = 1;           // Displacing solution concentration [mmol/kgw]
24
25  // Get 1D mimetic operators
26  Gradient G(k, m, dx);
27  Divergence D(k, m, dx);
28  Interpol I(m, 0.5);
29
30  // Allocate fields
31  vec C(m + 2); // Scalar field (concentrations)
32  vec V(m + 1); // Vector field (velocities)
33
34  // Impose initial conditions
35  C(0) = C0;
36  V.fill(vel);
37
38  // Hydrodynamic dispersion coefficient [m^2/year]
39  dis *= vel; // 75
40
41  // dt = dt/R (retardation)
42  dt /= R;
43
44  // Time integration loop
45  for (int i = 0; i <= iter; i++) {
46
47    // First-order forward-time scheme
48    C += dt * (D * (dis * (G * C)) - D * (V % (I * C)));
49
50    // Right boundary condition (reflection)
51    C(m + 1) = C(m);
52  }
53
54  // Spit out the new concentrations!
55  cout << C;
56
57  return 0;
58}

Elliptic Example (Laplacian)#

Elliptic 2D Example using Laplacian (examples/cpp/elliptic2D.cpp)#
 1/**
 2 * This example uses MOLE to solve a 2D BVP
 3 * It is the equivalent to examples_MATLAB/minimal_poisson2D.m
 4 * The output can be plotted via:
 5 * gnuplot> plot 'solution' matrix with image
 6 */
 7
 8#include "mole.h"
 9#include <iostream>
10
11int main() {
12  int k = 2; // Operators' order of accuracy
13  int m = 9; // Vertical resolution
14  int n = 9; // Horizontal resolution
15
16  // Get mimetic operators
17  Laplacian L(k, m, n, 1, 1);
18  RobinBC BC(k, m, 1, n, 1, 1, 0); // Dirichlet BC
19  L = L + BC;
20
21  // Build RHS for system of equations
22  mat rhs(m + 2, n + 2, fill::zeros);
23  rhs.row(0) = 100 * ones(1, n + 2); // Known value at the bottom boundary
24
25  // Solve the system of linear equations
26#ifdef EIGEN
27  // Use Eigen only if SuperLU (faster) is not available
28  vec sol = Utils::spsolve_eigen(L, vectorise(rhs));
29#else
30  vec sol = spsolve(L, vectorise(rhs)); // Will use SuperLU
31#endif
32
33  // Print out the solution
34  cout << reshape(sol, m + 2, n + 2);
35
36  return 0;
37}

Schrödinger Example (Complex Operators)#

Schrödinger 1D Example (examples/cpp/schrodinger1D.cpp)#
 1/**
 2 * This example uses MOLE to solve the 1D Schrodinger equation
 3 */
 4
 5#include "mole.h"
 6#include <algorithm>
 7#include <iostream>
 8
 9int main() {
10
11  int k = 4;   // Operators' order of accuracy
12  Real a = -5; // Left boundary
13  Real b = 5;  // Right boundary
14  int m = 500; // Number of cells
15  vec grid = linspace(a, b, m);
16  Real dx = grid(1) - grid(0); // Step size
17
18  // Get mimetic Laplacian operator
19  Laplacian L(k, m - 2, dx);
20
21  std::transform(grid.begin(), grid.end(), grid.begin(),
22                 [](Real x) { return x * x; });
23
24  sp_mat V(m, m); // Potential energy operator
25  V.diag(0) = grid;
26
27  // Hamiltonian
28  sp_mat H = -0.5 * (sp_mat)L + V;
29
30  cx_vec eigval;
31  eig_gen(eigval, (mat)H); // Compute eigenvalues
32
33  eigval = sort(eigval);
34
35  cout << "Energy levels = [ ";
36  for (int i = 0; i < 4; ++i)
37    cout << real(eigval(i) / eigval(0)) << ' ';
38  cout << "]\n";
39
40  return 0;
41}

Notes and Considerations#

Contents