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.

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#

Contents