field.cpp

Go to the documentation of this file.

/// @file field.cpp
/// @brief ElectricSolver and joule_heating implementations.
#include "field.hpp"
#include <algorithm>
#include <unordered_map>
 
namespace physics {
 
// ============================================================
// ElectricSolver
// ============================================================
 
num::SolverResult ElectricSolver::solve_potential(ScalarField3D& phi,
                                                   const ScalarField3D& sigma,
                                                   const std::vector<ElectrodeBC>& bcs,
                                                   double tol, int max_iter) {
    const num::Grid3D& g  = phi.grid();
    const num::Grid3D& sg = sigma.grid();
    const int nx = g.nx(), ny = g.ny(), nz = g.nz();
    const int N  = nx * ny * nz;
    const double inv_dx2 = 1.0 / (g.dx() * g.dx());
 
    auto flat = [&](int i, int j, int k) -> int {
        return k * ny * nx + j * nx + i;
    };
 
    std::unordered_map<int, double> bc_map;
    bc_map.reserve(bcs.size());
    for (const auto& e : bcs)
        bc_map[e.flat_idx] = static_cast<double>(e.voltage);
 
    constexpr int DI[6] = { 1,-1, 0, 0, 0, 0};
    constexpr int DJ[6] = { 0, 0, 1,-1, 0, 0};
    constexpr int DK[6] = { 0, 0, 0, 0, 1,-1};
 
    constexpr double penalty = 1e10;
    num::Vector b(N, 0.0);
 
    for (const auto& e : bcs) {
        b[e.flat_idx] = penalty * static_cast<double>(e.voltage);
        const int ei = e.flat_idx % nx;
        const int ej = (e.flat_idx / nx) % ny;
        const int ek = e.flat_idx / (nx * ny);
        for (int d = 0; d < 6; ++d) {
            int ni = ei + DI[d], nj = ej + DJ[d], nk = ek + DK[d];
            if (ni < 0 || ni >= nx || nj < 0 || nj >= ny || nk < 0 || nk >= nz)
                continue;
            int nidx = flat(ni, nj, nk);
            if (bc_map.count(nidx)) continue;
            double sigma_face = 0.5 * (sg(ei,ej,ek) + sg(ni,nj,nk));
            b[nidx] += sigma_face * inv_dx2 * static_cast<double>(e.voltage);
        }
    }
 
    auto matvec = [&](const num::Vector& v, num::Vector& Av) {
        for (int k = 0; k < nz; ++k)
        for (int j = 0; j < ny; ++j)
        for (int i = 0; i < nx; ++i) {
            int id = flat(i, j, k);
            if (bc_map.count(id)) { Av[id] = penalty * v[id]; continue; }
            double Av_ijk = 0.0;
            for (int d = 0; d < 6; ++d) {
                int ni = std::max(0, std::min(i + DI[d], nx-1));
                int nj = std::max(0, std::min(j + DJ[d], ny-1));
                int nk = std::max(0, std::min(k + DK[d], nz-1));
                int nidx = flat(ni, nj, nk);
                double sigma_face = 0.5 * (sg(i,j,k) + sg(ni,nj,nk));
                double v_nb = bc_map.count(nidx) ? 0.0 : v[nidx];
                Av_ijk += sigma_face * inv_dx2 * (v[id] - v_nb);
            }
            Av[id] = Av_ijk;
        }
    };
 
    num::Vector xv = phi.grid().to_vector();
    auto result = num::cg_matfree(matvec, b, xv, tol,
                                  static_cast<num::idx>(max_iter));
    phi.grid().from_vector(xv);
    return result;
}
 
ScalarField3D ElectricSolver::joule_heating(const ScalarField3D& sigma,
                                             const ScalarField3D& phi) {
    const num::Grid3D& g  = phi.grid();
    const num::Grid3D& sg = sigma.grid();
    const int nx = g.nx(), ny = g.ny(), nz = g.nz();
    const double inv2dx = 1.0 / (2.0 * g.dx());
 
    ScalarField3D Q(nx, ny, nz, phi.dx(), phi.ox(), phi.oy(), phi.oz());
    for (int k = 0; k < nz; ++k)
    for (int j = 0; j < ny; ++j)
    for (int i = 0; i < nx; ++i) {
        int ip = std::min(i+1,nx-1), im = std::max(i-1,0);
        int jp = std::min(j+1,ny-1), jm = std::max(j-1,0);
        int kp = std::min(k+1,nz-1), km = std::max(k-1,0);
        double dpdx = (g(ip,j,k) - g(im,j,k)) * inv2dx;
        double dpdy = (g(i,jp,k) - g(i,jm,k)) * inv2dx;
        double dpdz = (g(i,j,kp) - g(i,j,km)) * inv2dx;
        Q.grid().set(i, j, k, sg(i,j,k) * (dpdx*dpdx + dpdy*dpdy + dpdz*dpdz));
    }
    return Q;
}
 
} // namespace physics