probabilitycurrent2d_tests.cpp

Go to the documentation of this file.

/**
 * @file probabilitycurrent2d_tests.cpp
 * @brief Implements regression and invariant tests for probabilitycurrent2d tests.
 * @ingroup qpe_tests
 *
 * Source: @c tests/probabilitycurrent2d_tests.cpp. This file is part of Quantum Phenomena Explorer (QPE),
 * maintained by Officina Turini. See LICENSE_IT.txt / LICENSE.txt for the
 * QPE freeware source-available terms and THIRD_PARTY_NOTICES.txt for
 * third-party components.
 */
 
#include "probabilitycurrent2d_tests.h"
 
#include "core/quantum/probabilitycurrent2d.h"
 
#include <cmath>
#include <exception>
 
namespace qpe {
namespace tests {
 
/**
 * @brief Runs probability current2 d tests.
 * @param error Output diagnostic text on failure.
 */
bool runProbabilityCurrent2DTests(std::string &error)
{
    error.clear();
 
    try {
        // A real constant field has zero phase gradient and therefore J=0.
        ComplexField2D constantField(32, 24, 0.1f, 0.1f);
        constantField.fill(ComplexField2D::Value(1.0f, 0.0f));
        const ProbabilityCurrentField2D zeroCurrent =
            ProbabilityCurrent2D::compute(constantField, 1.0f, 1.0f);
        if (zeroCurrent.maxMagnitude > 1.0e-6f ||
            std::abs(zeroCurrent.integratedX) > 1.0e-8 ||
            std::abs(zeroCurrent.integratedY) > 1.0e-8) {
            error = "real constant field must have zero probability current";
            return false;
        }
 
        // Plane wave exp(i k x): the centred finite difference gives
        // Jx = |psi|^2 (hbar/m) sin(k dx)/dx in the interior.
        const std::size_t width = 96;
        const std::size_t height = 48;
        const float dx = 0.05f;
        const float dy = 0.05f;
        const double k = 1.4;
        ComplexField2D plane(width, height, dx, dy);
        for (std::size_t y = 0; y < height; ++y) {
            for (std::size_t x = 0; x < width; ++x) {
                const double phase = k * static_cast<double>(x) * dx;
                plane(x, y) = ComplexField2D::Value(
                    static_cast<float>(std::cos(phase)),
                    static_cast<float>(std::sin(phase)));
            }
        }
        plane.normalize();
 
        const float hbar = 1.0f;
        const float mass = 2.0f;
        const ProbabilityCurrentField2D current =
            ProbabilityCurrent2D::compute(plane, hbar, mass);
 
        const std::size_t cx = width / 2;
        const std::size_t cy = height / 2;
        const std::size_t ci = cy * width + cx;
        const double probability = std::norm(std::complex<double>(
            static_cast<double>(plane(cx, cy).real()),
            static_cast<double>(plane(cx, cy).imag())));
        const double expected = probability *
            (static_cast<double>(hbar) / static_cast<double>(mass)) *
            std::sin(k * dx) / dx;
 
        if (std::abs(static_cast<double>(current.jx[ci]) - expected) > 2.0e-5) {
            error = "plane-wave Jx does not match the finite-difference reference";
            return false;
        }
        if (std::abs(current.jy[ci]) > 1.0e-6f) {
            error = "plane wave along X must have Jy approximately zero";
            return false;
        }
        if (!(current.maxMagnitude > 0.0f) || !(current.integratedX > 0.0)) {
            error = "positive-k plane wave must carry positive integrated X current";
            return false;
        }
 
        // The gauge-covariant expression must reduce to the ordinary current
        // when every Peierls link is the identity.
        U1GaugeField2D identityGauge(width, height, dx, dy);
        const ProbabilityCurrentField2D gaugeCurrent =
            ProbabilityCurrent2D::computeGaugeCovariant(plane, identityGauge, hbar, mass);
        if (std::abs(static_cast<double>(gaugeCurrent.jx[ci]) -
                     static_cast<double>(current.jx[ci])) > 2.0e-7 ||
            std::abs(static_cast<double>(gaugeCurrent.jy[ci]) -
                     static_cast<double>(current.jy[ci])) > 2.0e-7 ||
            std::abs(gaugeCurrent.integratedX - current.integratedX) > 2.0e-7 ||
            std::abs(gaugeCurrent.integratedY - current.integratedY) > 2.0e-7) {
            error = "identity-gauge current must match the ordinary current";
            return false;
        }
 
 
        // Local U(1) gauge transformation: psi -> exp(i chi) psi and
        // U_ij -> exp(i chi_i) U_ij exp(-i chi_j).  The physical current
        // must remain unchanged.
        ComplexField2D transformed = plane;
        U1GaugeField2D transformedGauge(width, height, dx, dy);
        const double ax = 0.17;
        const double ay = -0.11;
        for (std::size_t y = 0; y < height; ++y) {
            for (std::size_t x = 0; x < width; ++x) {
                const double chi = ax * static_cast<double>(x) + ay * static_cast<double>(y);
                const std::complex<float> phase(static_cast<float>(std::cos(chi)),
                                                static_cast<float>(std::sin(chi)));
                transformed(x, y) *= phase;
                if (x + 1u < width) {
                    const double dchi = chi - (ax * static_cast<double>(x + 1u) + ay * static_cast<double>(y));
                    transformedGauge.setLinkX(x, y, U1GaugeField2D::Link(
                        static_cast<float>(std::cos(dchi)), static_cast<float>(std::sin(dchi))));
                }
                if (y + 1u < height) {
                    const double dchi = chi - (ax * static_cast<double>(x) + ay * static_cast<double>(y + 1u));
                    transformedGauge.setLinkY(x, y, U1GaugeField2D::Link(
                        static_cast<float>(std::cos(dchi)), static_cast<float>(std::sin(dchi))));
                }
            }
        }
        const ProbabilityCurrentField2D invariantCurrent =
            ProbabilityCurrent2D::computeGaugeCovariant(transformed, transformedGauge, hbar, mass);
        if (std::abs(static_cast<double>(invariantCurrent.jx[ci]) -
                     static_cast<double>(current.jx[ci])) > 3.0e-6 ||
            std::abs(static_cast<double>(invariantCurrent.jy[ci]) -
                     static_cast<double>(current.jy[ci])) > 3.0e-6 ||
            std::abs(invariantCurrent.integratedX - current.integratedX) > 3.0e-6 ||
            std::abs(invariantCurrent.integratedY - current.integratedY) > 3.0e-6) {
            error = "gauge-covariant probability current is not U(1) invariant";
            return false;
        }
    } catch (const std::exception &e) {
        error = e.what();
        return false;
    }
 
    return true;
}
 
} // namespace tests
} // namespace qpe