Boundary Conditions

MOLE supports a variety of boundary conditions for solving PDEs with different physical constraints.

Mixed Boundary Conditions

The MixedBC class implements mixed boundary conditions in the MOLE library.

API Reference

class MixedBC : public sp_mat

Mimetic Mixed Boundary Condition operator.

Public Functions

MixedBC(u16 k, u32 m, Real dx, const std::string &left, const std::vector<Real> &coeffs_left, const std::string &right, const std::vector<Real> &coeffs_right)

1-D Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells

  • dx – Spacing between cells

  • left – Type of boundary condition at the left boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_left – Coefficients for the left boundary condition

  • right – Type of boundary condition at the right boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_right – Coefficients for the right boundary condition

MixedBC(u16 k, u32 m, Real dx, u32 n, Real dy, const std::string &left, const std::vector<Real> &coeffs_left, const std::string &right, const std::vector<Real> &coeffs_right, const std::string &bottom, const std::vector<Real> &coeffs_bottom, const std::string &top, const std::vector<Real> &coeffs_top)

2-D Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells along x-axis

  • dx – Spacing between cells along x-axis

  • n – Number of cells along y-axis

  • dy – Spacing between cells along y-axis

  • left – Type of boundary condition at the left boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_left – Coefficients for the left boundary condition

  • right – Type of boundary condition at the right boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_right – Coefficients for the right boundary condition

  • bottom – Type of boundary condition at the bottom boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_bottom – Coefficients for the bottom boundary condition

  • top – Type of boundary condition at the top boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_top – Coefficients for the top boundary condition

MixedBC(u16 k, u32 m, Real dx, u32 n, Real dy, u32 o, Real dz, const std::string &left, const std::vector<Real> &coeffs_left, const std::string &right, const std::vector<Real> &coeffs_right, const std::string &bottom, const std::vector<Real> &coeffs_bottom, const std::string &top, const std::vector<Real> &coeffs_top, const std::string &front, const std::vector<Real> &coeffs_front, const std::string &back, const std::vector<Real> &coeffs_back)

3-D Constructor

Parameters:

  • k – Order of accuracy

  • m – Number of cells along x-axis

  • dx – Spacing between cells along x-axis

  • n – Number of cells along y-axis

  • dy – Spacing between cells along y-axis

  • o – Number of cells along z-axis

  • dz – Spacing between cells along z-axis

  • left – Type of boundary condition at the left boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_left – Coefficients for the left boundary condition

  • right – Type of boundary condition at the right boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_right – Coefficients for the right boundary condition

  • bottom – Type of boundary condition at the bottom boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_bottom – Coefficients for the bottom boundary condition

  • top – Type of boundary condition at the top boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_top – Coefficients for the top boundary condition

  • front – Type of boundary condition at the front boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_front – Coefficients for the front boundary condition

  • back – Type of boundary condition at the back boundary (‘Dirichlet’, ‘Neumann’, ‘Robin’)

  • coeffs_back – Coefficients for the back boundary condition

Robin Boundary Conditions

The RobinBC class implements Robin boundary conditions in the MOLE library.

API Reference

class RobinBC : public sp_mat

Mimetic Robin Boundary Condition operator.

Public Functions

RobinBC(u16 k, u32 m, Real dx, Real a, Real b)

1-D Robin boundary constructor

Parameters:

  • k – mimetic order of accuracy

  • m – number of cells in x-dimension

  • dx – cell width in x-direction

  • a – Coefficient of the Dirichlet function

  • b – Coefficient of the Neumann function

RobinBC(u16 k, u32 m, Real dx, u32 n, Real dy, Real a, Real b)

2-D Robin boundary constructor

Note

Uses 1-D Robin to build the 2-D operator

Parameters:

  • k – mimetic order of accuracy

  • m – number of cells in x-dimension

  • dx – cell width in x-direction

  • n – number of cells in y-dimension

  • dy – cell width in y-direction

  • a – Coefficient of the Dirichlet function

  • b – Coefficient of the Neumann function

RobinBC(u16 k, u32 m, Real dx, u32 n, Real dy, u32 o, Real dz, Real a, Real b)

3-D Robin boundary constructor

Note

Uses 1-D Robin to build the 3-D operator

Parameters:

  • k – mimetic order of accuracy

  • m – number of cells in x-dimension

  • dx – cell width in x-direction

  • n – number of cells in y-dimension

  • dy – cell width in y-direction

  • o – number of cells in z-dimension

  • dz – cell width in z-direction

  • a – Coefficient of the Dirichlet function

  • b – Coefficient of the Neumann function

Usage Examples

Elliptic Problem with Mixed Boundary Conditions

Elliptic 1D Example with Mixed Boundary Conditions (examples/cpp/elliptic1D.cpp)
 1/**
 2 * This example uses MOLE to solve a 1D BVP
 3 */
 4
 5#include "mole.h"
 6#include <iostream>
 7
 8int main() {
 9
10  int k = 6;             // Operators' order of accuracy
11  Real a = 0;            // Left boundary
12  Real b = 1;            // Right boundary
13  int m = 2 * k + 1;     // Number of cells
14  Real dx = (b - a) / m; // Step size
15
16  // Get mimetic operators
17  Laplacian L(k, m, dx);
18  Real d = 1; // Dirichlet coefficient
19  Real n = 1; // Neumann coefficient
20  RobinBC BC(k, m, dx, d, n);
21  L = L + BC;
22
23  // 1D Staggered grid
24  vec grid(m + 2);
25  grid(0) = a;
26  grid(1) = a + dx / 2.0;
27  int i;
28  for (i = 2; i <= m; i++) {
29    grid(i) = grid(i - 1) + dx;
30  }
31  grid(i) = b;
32
33  // Build RHS for system of equations
34  vec rhs(m + 2);
35  rhs = exp(grid); // rhs(0) = 1
36  rhs(0) = 0;
37  rhs(m + 1) = 2 * exp(1); // rhs(1) = 2e
38
39  // Solve the system of linear equations
40#ifdef EIGEN
41  // Use Eigen only if SuperLU (faster) is not available
42  vec sol = Utils::spsolve_eigen(L, rhs);
43#else
44  vec sol = spsolve(L, rhs); // Will use SuperLU
45#endif
46
47  // Print out the solution
48  cout << sol;
49
50  return 0;
51}

Complex Problem with Advanced Boundary Handling

3D Convection-Diffusion Example with Complex Boundary Conditions (examples/cpp/convection_diffusion3D.cpp)
  1/**
  2 * @file convection_diffusion3D.cpp
  3 * @brief Solves a 3D Convection-Diffusion equation using Mimetic Operators for Linear Equations (MOLE).
  4 *
  5 * This program numerically solves the 3D convection-diffusion equation using MOLE techniques.
  6 * It simulates the distribution of CO2 concentration over time within a porous medium,
  7 * considering both diffusion and advection effects. 
  8 * 
  9 * If OUTPUT_FRAME_DATA is set to 1, slices (2D cross-sections) of the concentration field are extracted 
 10 * and saved to a text file for visualization.
 11 */
 12
 13#include <iostream>
 14#include <fstream>
 15#include <mole.h>
 16
 17#define OUTPUT_FRAME_DATA 0
 18
 19// Helper for linear indexing in 3D arrays (0-based indexing)
 20inline size_t idx3D(size_t i, size_t j, size_t k, size_t m, size_t n, size_t o) {
 21    return i + (m+2)*j + (m+2)*(n+2)*k;
 22}
 23
 24int main() {
 25    // Parameters
 26    unsigned short k = 2;
 27    unsigned int m = 101;
 28    unsigned int n = 51;
 29    unsigned int o = 101;
 30
 31    double a = 0.0, b = 101.0;
 32    double c = 0.0, d = 51.0;
 33    double e = 0.0, f = 101.0;
 34
 35    double dx = (b - a) / m;
 36    double dy = (d - c) / n;
 37    double dz = (f - e) / o;
 38
 39    // Construct operators
 40    sp_mat D = Divergence(k, m, n, o, dx, dy, dz);
 41    sp_mat G = Gradient(k, m, n, o, dx, dy, dz);
 42    sp_mat I = Interpol(m, n, o, 1, 1, 1);
 43
 44    size_t scalarSize = (m+2)*(n+2)*(o+2);
 45    size_t vectorSize = G.n_rows;
 46
 47    // Allocate fields
 48    std::vector<double> V(vectorSize, 0.0);
 49    vec C(scalarSize, fill::zeros);
 50
 51    // Initial conditions
 52    int bottom = 10; 
 53    int top = 15;   
 54    int seal = 40;  
 55    int seal_idx = seal - 1;
 56    int seal5_idx = (seal + 5) - 1;
 57
 58    // Construct the velocity field
 59    size_t yCount = m*(n+1)*o;
 60    std::vector<double> y(yCount, 1.0);
 61
 62    // Apply shale conditions on velocity field
 63    for (int i_ = 0; i_ < (int)m; i_++) {
 64        for (int k_ = 0; k_ < (int)o; k_++) {
 65            y[i_ + m*seal_idx + m*(n+1)*k_] = 0.0;
 66            y[i_ + m*seal5_idx + m*(n+1)*k_] = 0.0;
 67        }
 68    }
 69
 70    // Assign y into V at the correct offset
 71    size_t offset = (m+1)*n*o; 
 72    for (size_t i_ = 0; i_ < yCount; ++i_) {
 73        if (offset + i_ < V.size()) {
 74            V[offset + i_] = y[i_];
 75        }
 76    }
 77
 78    // Set initial density
 79    int mid_x = (int)std::ceil((m+2)/2.0) - 1; 
 80    int mid_z = (int)std::ceil((o+2)/2.0) - 1;
 81    for (int j = bottom - 1; j <= top - 1; j++) {
 82        size_t idx = idx3D(mid_x, j, mid_z, m, n, o);
 83        C(idx) = 1.0;
 84    }
 85
 86    // Well indices where C=1
 87    std::vector<size_t> wellIndices;
 88    for (size_t i_ = 0; i_ < scalarSize; ++i_) {
 89        if (C(i_) == 1.0) {
 90            wellIndices.push_back(i_);
 91        }
 92    }
 93
 94    // Diffusivity and porosity
 95    double diff = 1.0;
 96    double porosity = 1.0;
 97    diff *= porosity;
 98
 99    // Build K
100    std::vector<double> K(vectorSize, diff);
101    std::vector<double> kk(yCount, diff);
102    for (int i_ = 0; i_ < (int)m; i_++) {
103        for (int k_ = 0; k_ < (int)o; k_++) {
104            kk[i_ + m*seal_idx + m*(n+1)*k_] = diff/10.0;
105            kk[i_ + m*seal5_idx + m*(n+1)*k_] = diff/40.0;
106        }
107    }
108    for (size_t i_ = 0; i_ < yCount; ++i_) {
109        if (offset + i_ < K.size()) {
110            K[offset + i_] = kk[i_];
111        }
112    }
113
114    // Time step calculation
115    double dt1 = dx*dx/(3*diff)/3.0;
116    double maxV = 0.0;
117    for (auto val : V) {
118        if (val > maxV) maxV = val;
119    }
120    double dt2 = (maxV > 0.0) ? (dx/maxV)/3.0 : 1e-3;
121    double dt = std::min(dt1, dt2);
122    int iters = 120;
123
124    // Convert V and K to arma::vec
125    arma::vec K_arma(K);
126    arma::vec V_arma(V);
127
128    arma::ivec offsets_vec(1);
129    offsets_vec(0) = 0;
130
131    sp_mat Kdiag = spdiags(K_arma, offsets_vec, K_arma.n_elem, K_arma.n_elem);
132    sp_mat Vdiag = spdiags(V_arma, offsets_vec, V_arma.n_elem, V_arma.n_elem);
133
134    SizeMat size_identity(D.n_rows, D.n_rows);
135    sp_mat I_sp = speye(size_identity);
136
137    // Operators: L and Dadv
138    sp_mat L = dt * D * Kdiag * G + I_sp;
139    sp_mat Dadv = dt * D * Vdiag * I;
140
141    #if OUTPUT_FRAME_DATA
142    // Open a single file to store selected frames
143    std::ofstream frameFile("frames.txt");
144    if(!frameFile) {
145        std::cerr << "Error opening frames.txt for writing.\n";
146        return 1;
147    }
148    #endif
149
150    // Time-stepping loop
151    for (int i_ = 1; i_ <= iters*3; ++i_) {
152        // Diffusion step
153        vec Cnew = L * C;
154        for (auto w : wellIndices) {
155            Cnew(w) = 1.0;
156        }
157        C = Cnew;
158
159        // Advection step
160        vec Cadv = Dadv * C;
161        Cadv = C - Cadv;
162        for (auto w : wellIndices) {
163            Cadv(w) = 1.0;
164        }
165        C = Cadv;
166
167        #if OUTPUT_FRAME_DATA
168        // Write only selected frames to a single file
169        frameFile << "FRAME " << i_ << "\n";
170        for (int j = 0; j < (int)(n+2); j++) {
171            for (int k_ = 0; k_ < (int)(o+2); k_++) {
172                size_t idx = idx3D(seal_idx, j, k_, m, n, o);
173                frameFile << C(idx);
174                if (k_ < (int)(o+1)) frameFile << " ";
175            }
176            frameFile << "\n";
177        }
178        frameFile << "\n"; // Blank line between frames
179        #endif
180    }
181    
182    #if OUTPUT_FRAME_DATA
183    frameFile.close();
184    #endif
185
186    // Display minimum and maximum concentration values
187    std::cout << "Minimum CO2 concentration: " << C.min() << "\n";
188    std::cout << "Maximum CO2 concentration: " << C.max() << "\n";
189
190    return 0;
191}

Importance of Boundary Conditions

Boundary conditions are critical for ensuring that differential equation solutions are unique and physically meaningful. They specify constraints at the boundaries of the computational domain.