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.

Public Functions

Gradient(u16 k, u32 m, Real dx)

1-D Mimetic Gradient Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells

  • dx – Spacing between cells

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

2-D Mimetic Gradient 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

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

3-D Mimetic Gradient 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

vec getP()

Returns the weights used in the Mimeitc Gradient Operators.

Note

for informational purposes only, can be used in constructing new operators.

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.

Public Functions

Divergence(u16 k, u32 m, Real dx)

1-D Mimetic Divergence Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells

  • dx – Spacing between cells

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

2-D Mimetic Divergence 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

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

3-D Mimetic Divergence 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

vec getQ()

Returns the weights used in the Mimeitc Divergence Operators.

Note

for informational purposes only, can be used in constructing new operators.

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