schrod_eig.cpp

Go to the documentation of this file.

/**
 * \file schrod_eig.cpp
 * \example schrod _eig.cpp
 *
 * Calculate wave functions of Schrödinger equation in 2d problems.
 *
 */
#ifndef EXECUTABLE_DIMENSION
  #define EXECUTABLE_DIMENSION 2
#endif
 
#include <MoFEM.hpp>
#include <schrod_eig.hpp>
#undef EPS
#include <slepceps.h>
 
using namespace MoFEM;
using namespace schrod_eig;
 
template <int DIM> struct ElementsAndOps {};
 
template <> struct ElementsAndOps<2> {
  using DomainEle = FaceElementForcesAndSourcesCore;
};
 
template <> struct ElementsAndOps<3> {
  using DomainEle = VolumeElementForcesAndSourcesCore;
};
 
constexpr int SPACE_DIM =
    EXECUTABLE_DIMENSION; // Space dimension of problem, mesh
 
using DomainEle = ElementsAndOps<SPACE_DIM>::DomainEle;
using DomainEleOp = DomainEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
 
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpGradGrad<1, 1, SPACE_DIM>;
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpMass<1, 1>;
 
//! [Physical constants and parameters]
double hbar = 1.054571817e-34; // Reduced Planck Constant [J*s]
double m0 = 9.1093837015e-31;  // Free electron mass [kg]
double scale = 1e-9;           // Length scaling factor: nanometer → meter [m]
double effMass = 1;            // Effective mass in units of m0 [a.u.]
double q = 1.602176634e-19;    // Elementary charge [C]
double potential = 0 * q;      // Potential energy inside the potential well [J]
//! [Physical constants and parameters]
 
int order = 2;
 
static char help[] = "...\n\n";
 
struct Example {
 
  Example(MoFEM::Interface &m_field) : mField(m_field) {}
 
  MoFEMErrorCode runProblem();
 
private:
  MoFEM::Interface &mField;
  Simple *simple;
 
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode outputResults();
  MoFEMErrorCode checkResults();
 
  SmartPetscObj<Mat> M; // Mass matrix
  SmartPetscObj<Mat> H; // Hamiltonian matrix
  SmartPetscObj<EPS> eps;
};
 
//! [Run problem]
MoFEMErrorCode Example::runProblem() {
  MoFEMFunctionBegin;
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  CHKERR outputResults();
  CHKERR checkResults();
  MoFEMFunctionReturn(0);
}
//! [Run problem]
 
//! [Read mesh]
MoFEMErrorCode Example::readMesh() {
  MoFEMFunctionBegin;
  simple = mField.getInterface<Simple>();
 
  MOFEM_LOG("EXAMPLE", Sev::inform)
      << "Read mesh for problem in " << EXECUTABLE_DIMENSION;
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  MoFEMFunctionReturn(0);
}
//! [Read mesh]
 
//! [Set up problem]
MoFEMErrorCode Example::setupProblem() {
  MoFEMFunctionBegin;
  // Add field
  CHKERR simple->addDomainField("PHI", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE, 1);
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
  CHKERR simple->setFieldOrder("PHI", order);
  CHKERR simple->setUp();
  MoFEMFunctionReturn(0);
}
//! [Set up problem]
 
//! [Applying essential BC]
MoFEMErrorCode Example::boundaryCondition() {
  MoFEMFunctionBegin;
 
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
      simple->getProblemName(), "BOUNDARY", std::string("PHI"),
      true); // Dirichlet BC. (PHI=0 on the boundaries).
 
  MoFEMFunctionReturn(0);
}
//! [Applying essential BC]
 
//! [Push operators to pipeline]
MoFEMErrorCode Example::assembleSystem() {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-effMass", &effMass,
                               PETSC_NULLPTR);
 
  auto get_kinetic_coe = [&](const double, const double, const double) {
    return hbar * hbar /
           (2.0 * effMass * m0 *
            (scale * scale)); // Coefficient for kinetic energy term
  };
  auto get_potential = [](const double, const double, const double) {
    return potential;
  };
  auto get_mass_coe = [](const double, const double, const double) {
    return 1.0;
  };
 
  auto dm = simple->getDM();
  CHKERR DMCreateMatrix_MoFEM(dm, H);
  M = matDuplicate(H, MAT_SHARE_NONZERO_PATTERN);
 
  auto calculate_Hamiltonian = [&]() {
    MoFEMFunctionBegin;
    pipeline_mng->getDomainLhsFE().reset();
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pipeline_mng->getOpDomainLhsPipeline(), {H1});
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpDomainGradGrad("PHI", "PHI", get_kinetic_coe));
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpDomainMass("PHI", "PHI", get_potential));
    auto integration_rule = [](int, int, int approx_order) {
      return 2 * (approx_order - 1);
    };
 
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
    pipeline_mng->getDomainLhsFE()->B = H;
    CHKERR MatZeroEntries(H);
    CHKERR pipeline_mng->loopFiniteElements();
    CHKERR MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_mass_matrix = [&]() {
    MoFEMFunctionBegin;
    pipeline_mng->getDomainLhsFE().reset();
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pipeline_mng->getOpDomainLhsPipeline(), {H1});
    pipeline_mng->getOpDomainLhsPipeline().push_back(
        new OpDomainMass("PHI", "PHI", get_mass_coe));
 
    auto integration_rule = [](int, int, int approx_order) {
      return 2 * approx_order;
    };
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
    CHKERR MatZeroEntries(M);
    pipeline_mng->getDomainLhsFE()->B = M;
    CHKERR pipeline_mng->loopFiniteElements();
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR calculate_Hamiltonian();
  CHKERR calculate_mass_matrix();
 
  MoFEMFunctionReturn(0);
}
//! [Push operators to pipeline]
 
//! [Solve]
MoFEMErrorCode Example::solveSystem() {
  MoFEMFunctionBegin;
 
  auto create_eps = [](MPI_Comm comm) {
    EPS eps;
    CHKERR EPSCreate(comm, &eps);
    return SmartPetscObj<EPS>(eps);
  };
 
  auto setup_eps = [&]() {
    MoFEMFunctionBegin;
    CHKERR EPSSetProblemType(eps, EPS_GHEP);
    CHKERR EPSSetWhichEigenpairs(eps, EPS_SMALLEST_MAGNITUDE);
    CHKERR EPSSetFromOptions(eps);
    PetscInt nev = 20;
    EPSSetDimensions(eps, nev, PETSC_DEFAULT, PETSC_DEFAULT);
    MoFEMFunctionReturn(0);
  };
 
  auto print_info = [&]() {
    MoFEMFunctionBegin;
    ST st;
    EPSType type;
    PetscReal tol;
    PetscInt nev, maxit, its;
    // Optional: Get some information from the solver and display it
    CHKERR EPSGetIterationNumber(eps, &its);
    MOFEM_LOG_C("EXAMPLE", Sev::inform,
                " Number of iterations of the method: %d", its);
    CHKERR EPSGetST(eps, &st);
    CHKERR EPSGetType(eps, &type);
    MOFEM_LOG_C("EXAMPLE", Sev::inform, " Solution method: %s", type);
    CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
    MOFEM_LOG_C("EXAMPLE", Sev::inform, " Number of requested eigenvalues: %d",
                nev);
    CHKERR EPSGetTolerances(eps, &tol, &maxit);
    MOFEM_LOG_C("EXAMPLE", Sev::inform,
                " Stopping condition: tol=%.4g, maxit=%d", (double)tol, maxit);
 
    MoFEMFunctionReturn(0);
  };
 
  // Create eigensolver context
  eps = create_eps(mField.get_comm());
  CHKERR EPSSetOperators(eps, H, M);
 
  // Setup EPS
  CHKERR setup_eps();
 
  // Solve problem
  CHKERR EPSSolve(eps);
 
  // Print info
  CHKERR print_info();
 
  MoFEMFunctionReturn(0);
}
//! [Solve]
 
//! [Postprocess results]
MoFEMErrorCode Example::outputResults() {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  pipeline_mng->getDomainLhsFE().reset();
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
  auto phi_ptr = boost::make_shared<VectorDouble>();
  auto square_ptr = boost::make_shared<VectorDouble>();
 
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues("PHI", phi_ptr));
  post_proc_fe->getOpPtrVector().push_back(new OpSquare(
      square_ptr,
      phi_ptr)); // Calculate square of wave function for checking probability density.
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(post_proc_fe->getPostProcMesh(),
                  post_proc_fe->getMapGaussPts(),
 
                  OpPPMap::DataMapVec{{"PHI", phi_ptr}, {"SQUARE", square_ptr}},
 
                  OpPPMap::DataMapMat{},
 
                  OpPPMap::DataMapMat{},
 
                  OpPPMap::DataMapMat{})
 
  );
 
  pipeline_mng->getDomainRhsFE() = post_proc_fe;
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
 
  PetscInt nev, nconv, n_output = 0;
  CHKERR EPSGetDimensions(eps, &nev, PETSC_NULLPTR, PETSC_NULLPTR);
  CHKERR EPSGetConverged(eps, &nconv);
  n_output = std::min(nconv, nev);
  if (nconv < nev) {
    MOFEM_LOG_C("EXAMPLE", Sev::warning,
                " Only %" PetscInt_FMT " of %" PetscInt_FMT
                " requested eigenpairs converged",
                nconv, nev);
  }
  PetscScalar eigr, eigi;
  for (PetscInt nn = 0; nn < n_output; nn++) {
    CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, D, PETSC_NULLPTR);
    CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
    MOFEM_LOG_C("EXAMPLE", Sev::inform,
                " Eigenpair = %" PetscInt_FMT " Eigen Energy = %.8g eV", nn,
                eigr / q); // Convert the unit Joule to eV for output.
    CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
    CHKERR pipeline_mng->loopFiniteElements();
    post_proc_fe->writeFile("out_schrod_" +
                            boost::lexical_cast<std::string>(nn) + ".h5m");
  }
 
  MoFEMFunctionReturn(0);
}
//! [Postprocess results]
 
//! [Check]
MoFEMErrorCode Example::checkResults() {
  MoFEMFunctionBegin;
  PetscBool test_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test", &test_flg,
                             PETSC_NULLPTR);
  if (test_flg) {
    PetscInt nconv;
    CHKERR EPSGetConverged(eps, &nconv);
    if (nconv < 1) {
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "No eigenpairs converged");
    }
    PetscScalar eigr, eigi;
    CHKERR EPSGetEigenpair(eps, 0, &eigr, &eigi, PETSC_NULLPTR, PETSC_NULLPTR);
    constexpr double regression_value =
        0.188; // Check the result for ground energy level [eV]
    if (fabs(eigr / q - regression_value) > 1e-3) {
      PetscPrintf(PETSC_COMM_WORLD,
                  "Calculated ground energy: %.4g, expected %.4g\n",
                  (double)eigr / q, regression_value);
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Regression test faileed; wrong eigen value. Try higher order or "
              "finer mesh.");
    }
  }
  MoFEMFunctionReturn(0);
}
//! [Check]
 
int main(int argc, char *argv[]) {
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  SlepcInitialize(&argc, &argv, param_file, help);
  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
 
  // Add logging channel for example
  auto core_log = logging::core::get();
  core_log->add_sink(
      LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
  LogManager::setLog("EXAMPLE");
  MOFEM_LOG_TAG("EXAMPLE", "example");
 
  try {
 
    //! [Register MoFEM discrete manager in PETSc]
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    //! [Register MoFEM discrete manager in PETSc
 
    //! [Create MoAB]
    moab::Core mb_instance;              ///< mesh database
    moab::Interface &moab = mb_instance; ///< mesh database interface
    //! [Create MoAB]
 
    //! [Create MoFEM]
    MoFEM::Core core(moab);           ///< finite element database
    MoFEM::Interface &m_field = core; ///< finite element database insterface
    //! [Create MoFEM]
 
    //! [Example]
    Example ex(m_field);
    CHKERR ex.runProblem();
    //! [Example]
  }
  CATCH_ERRORS;
 
  SlepcFinalize();
  CHKERR MoFEM::Core::Finalize();
}