initial_diffusion.cpp

Go to the documentation of this file.

/**
 * \file inital_diffusion.cpp
 * \example mofem/tutorials/scl-10/initial_diffusion.cpp
 *
 **/
 
#include <stdlib.h>
#include <cmath>
#include <MoFEM.hpp>
#include <SourceFunction.hpp>
 
#define BOOST_MATH_GAUSS_NO_COMPUTE_ON_DEMAND
 
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
template <int DIM> struct ElementsAndOps {};
 
//! [Define dimension]
constexpr int SPACE_DIM = 3; //< Space dimension of problem, mesh
//! [Define dimension]
 
using DomainEle = VolumeElementForcesAndSourcesCore;
using DomainEleOp = DomainEle::UserDataOperator;
using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
 
using OpDomainMass = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpMass<1, 1>;
using OpDomainGradGrad = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::BiLinearForm<GAUSS>::OpGradGrad<1, 1, SPACE_DIM>;
using OpDomainTimesScalarField = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::LinearForm<GAUSS>::OpBaseTimesScalar<1>;
using OpDomainGradTimesVec = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::LinearForm<GAUSS>::OpGradTimesTensor<1, 1, SPACE_DIM>;
using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
    PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
 
const double n = 1.44;  ///< refractive index of diffusive medium
const double c = 30.;   ///< speed of light (cm/ns)
const double v = c / n; ///< phase velocity of light in medium (cm/ns)
 
double mu_a;  ///< absorption coefficient (cm^-1)
double mu_sp; ///< scattering coefficient (cm^-1)
double D;
 
double slab_thickness;
double beam_radius; //< spot radius
double beam_centre_x;
double beam_centre_y;
double flux_magnitude = 1e3; ///< impulse magnitude
const int kronrod_points =
    15; ///< number of points for kronrod integration, can be 15, 31, 41, 51, or 61 (from boost library docs)
///< This has been tested and gives the same result for any number of points. Increasing the number of points will increase the compute time, so 15 is used as default.
double initial_time;
 
char out_file_name[255] = "init_file.dat";
;
int numHoLevels = 1;
 
PetscBool output_volume = PETSC_FALSE;
PetscBool testing = PETSC_FALSE;
 
#include <boost/math/quadrature/gauss_kronrod.hpp>
using namespace boost::math::quadrature;
 
struct PhotonDiffusion {
public:
  PhotonDiffusion(MoFEM::Interface &m_field);
 
  // Declaration of the main function to run analysis
  MoFEMErrorCode runProgram();
 
  static inline double sourceFunction(const double x, const double y,
                                      const double z) {
    return SourceFunction::sourceFunctionEval(
        x, y, z, beam_radius, beam_centre_x, beam_centre_y, slab_thickness,
        mu_a, mu_sp, flux_magnitude, initial_time, v, D);
  }
  //! [sourceFunction]
 
private:
  // Declaration of other main functions called in runProgram()
  MoFEMErrorCode readMesh();
  MoFEMErrorCode setupProblem();
  MoFEMErrorCode createCommonData();
  MoFEMErrorCode initialCondition();
  MoFEMErrorCode boundaryCondition();
  MoFEMErrorCode assembleSystem();
  MoFEMErrorCode solveSystem();
  MoFEMErrorCode checkResults();
  MoFEMErrorCode outputResults();
 
  // Main interfaces
  MoFEM::Interface &mField;
 
  struct CommonData {
    boost::shared_ptr<VectorDouble> uAtPtsPtr;
    SmartPetscObj<Vec> L2Vec;
    SmartPetscObj<Vec> resVec;
  };
 
  boost::shared_ptr<CommonData> commonDataPtr;
 
  struct OpError;
};
 
PhotonDiffusion::PhotonDiffusion(MoFEM::Interface &m_field) : mField(m_field) {}
 
struct PhotonDiffusion::OpError : public DomainEleOp {
  boost::shared_ptr<CommonData> commonDataPtr;
 
  OpError(boost::shared_ptr<CommonData> &common_data_ptr)
      : DomainEleOp("PHOTON_FLUENCE_RATE", OPROW),
        commonDataPtr(common_data_ptr) {}
 
  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
    MoFEMFunctionBegin;
 
    if (const size_t nb_dofs = data.getIndices().size()) {
 
      const int nb_integration_pts = getGaussPts().size2();
      auto t_w = getFTensor0IntegrationWeight();
      auto t_val = getFTensor0FromVec(*(commonDataPtr->uAtPtsPtr));
      auto t_coords = getFTensor1CoordsAtGaussPts();
 
      VectorDouble nf(nb_dofs, false);
      nf.clear();
 
      FTensor::Index<'i', 3> i;
      const double volume = getMeasure();
 
      auto t_row_base = data.getFTensor0N();
      double error = 0;
      for (int gg = 0; gg != nb_integration_pts; ++gg) {
 
        const double alpha = t_w * volume;
        double diff = t_val - PhotonDiffusion::sourceFunction(
                                  t_coords(0), t_coords(1), t_coords(2));
 
        error += alpha * pow(diff, 2);
 
        for (size_t r = 0; r != nb_dofs; ++r) {
          nf[r] += alpha * t_row_base * diff;
          ++t_row_base;
        }
 
        ++t_w;
        ++t_val;
        ++t_coords;
      }
 
      const int index = 0;
      CHKERR VecSetValue(commonDataPtr->L2Vec, index, error, ADD_VALUES);
      CHKERR VecSetValues(commonDataPtr->resVec, data, &nf[0], ADD_VALUES);
    }
 
    MoFEMFunctionReturn(0);
  }
};
 
MoFEMErrorCode PhotonDiffusion::readMesh() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR mField.getInterface(simple);
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::setupProblem() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
  CHKERR simple->addDomainField("PHOTON_FLUENCE_RATE", H1,
                                AINSWORTH_LEGENDRE_BASE, 1);
  CHKERR simple->addBoundaryField("PHOTON_FLUENCE_RATE", H1,
                                  AINSWORTH_LEGENDRE_BASE, 1);
 
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-flux_magnitude",
                               &flux_magnitude, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-slab_thickness",
                               &slab_thickness, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-beam_radius", &beam_radius,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-beam_centre_x",
                               &beam_centre_x, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-beam_centre_y",
                               &beam_centre_y, PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_a", &mu_a,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-mu_sp", &mu_sp,
                               PETSC_NULLPTR);
  CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-initial_time",
                               &initial_time, PETSC_NULLPTR);
 
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, "", "-output_file", out_file_name,
                               255, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-output_volume",
                             &output_volume, PETSC_NULLPTR);
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-testing", &testing,
                             PETSC_NULLPTR);
 
  D = 1. / (3. * (mu_a + mu_sp));
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Refractive index: " << n;
  MOFEM_LOG("INITIAL", Sev::inform) << "Speed of light (cm/ns): " << c;
  MOFEM_LOG("INITIAL", Sev::inform)
      << "Phase velocity in medium (cm/ns): " << v;
  MOFEM_LOG("INITIAL", Sev::inform)
      << "Absorption coefficient (cm^-1): " << mu_a;
  MOFEM_LOG("INITIAL", Sev::inform)
      << "Scattering coefficient (cm^-1): " << mu_sp;
  MOFEM_LOG("INITIAL", Sev::inform) << "Diffusion coefficient D : " << D;
  MOFEM_LOG("INITIAL", Sev::inform) << "Impulse magnitude: " << flux_magnitude;
  MOFEM_LOG("INITIAL", Sev::inform) << "Compute time (ns): " << initial_time;
  MOFEM_LOG("INITIAL", Sev::inform) << "Slab thickness: " << slab_thickness;
 
  int order = 2;
  CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Approximation order: " << order;
  MOFEM_LOG("INITIAL", Sev::inform) << "Kronrod points: " << kronrod_points;
 
  CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order);
 
  // if (numHoLevels > 0) {
 
  //   Range ho_ents;
  //   for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
  //     if (it->getName().compare(0, 3, "CAM") == 0) {
  //       CHKERR mField.get_moab().get_entities_by_dimension(it->getMeshset(), 2,
  //                                                          ho_ents, true);
  //     }
  //   }
 
  //   EntityHandle meshset;
  //   CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
  //   CHKERR mField.get_moab().add_entities(meshset, ho_ents);
  //   std::string field_name;
  //   field_name = "out_test_" +
  //                boost::lexical_cast<std::string>(mField.get_comm_rank()) +
  //                ".vtk";
  //   CHKERR mField.get_moab().write_file(field_name.c_str(), "VTK", "", &meshset,
  //                                       1);
  //   CHKERR mField.get_moab().delete_entities(&meshset, 1);
 
  //   CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(ho_ents);
 
  //   CHKERR simple->setFieldOrder("PHOTON_FLUENCE_RATE", order + 1, &ho_ents);
 
  //   CHKERR mField.getInterface<CommInterface>()->synchroniseFieldEntities(
  //       "PHOTON_FLUENCE_RATE");
  // }
 
  CHKERR simple->setUp();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::createCommonData() {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
 
  commonDataPtr = boost::make_shared<PhotonDiffusion::CommonData>();
  commonDataPtr->resVec = createDMVector(simple->getDM());
  commonDataPtr->L2Vec =
      createVectorMPI(mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
  commonDataPtr->uAtPtsPtr = boost::make_shared<VectorDouble>();
 
  MoFEMFunctionReturn(0);
}
 
MoFEMErrorCode PhotonDiffusion::initialCondition() {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}
 
//![interior_boundary]
MoFEMErrorCode PhotonDiffusion::boundaryCondition() {
  MoFEMFunctionBegin;
 
  auto *simple = mField.getInterface<Simple>();
 
  // Get boundary faces marked in block named "INT"
  Range boundary_faces;
  for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, BLOCKSET, it)) {
    std::string entity_name = it->getName();
    if (entity_name.compare(0, 3, "INT") == 0) {
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
                                                 boundary_faces, true);
    }
  }
 
  if (boundary_faces.empty()) {
    for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(mField, SIDESET, it)) {
      std::string entity_name = it->getName();
      CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
                                                 boundary_faces, true);
    }
  }
 
  // Get boundary edges in "INT"
  Range boundary_ents;
  CHKERR mField.get_moab().get_adjacencies(
      boundary_faces, 1, false, boundary_ents, moab::Interface::UNION);
  // Add vertices to boundary entities
  Range boundary_verts;
  CHKERR mField.get_moab().get_adjacencies(
      boundary_faces, 0, false, boundary_verts, moab::Interface::UNION);
 
  boundary_faces.merge(boundary_verts);
  boundary_faces.merge(boundary_ents);
 
  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
      boundary_faces);
 
  EntityHandle meshset;
  CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
  CHKERR mField.get_moab().add_entities(meshset, boundary_faces);
 
  // Remove DOFs as homogeneous boundary condition is used
 
  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
      simple->getProblemName(), "PHOTON_FLUENCE_RATE", boundary_faces);
 
  MoFEMFunctionReturn(0);
}
//! [interior_boundary]
 
//! [assembleSystem]
MoFEMErrorCode PhotonDiffusion::assembleSystem() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto integration_rule = [](int o_row, int o_col, int approx_order) {
    return 2 * approx_order;
  };
 
  auto set_domain = [&]() {
    MoFEMFunctionBegin;
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pipeline_mng->getOpDomainLhsPipeline(), {H1});
 
    pipeline_mng->getOpDomainLhsPipeline().push_back(new OpDomainMass(
        "PHOTON_FLUENCE_RATE", "PHOTON_FLUENCE_RATE",
        [](const double, const double, const double) { return 1; }));
 
    pipeline_mng->getOpDomainRhsPipeline().push_back(
        new OpDomainSource("PHOTON_FLUENCE_RATE", sourceFunction));
 
    CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
    MoFEMFunctionReturn(0);
  };
 
  auto set_boundary = [&]() {
    MoFEMFunctionBegin;
    CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
    CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
    MoFEMFunctionReturn(0);
  };
 
  CHKERR set_domain();
  CHKERR set_boundary();
 
  MoFEMFunctionReturn(0);
}
//! [assembleSystem]
 
//! [solveSystem]
MoFEMErrorCode PhotonDiffusion::solveSystem() {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto solver = pipeline_mng->createKSP();
 
  CHKERR KSPSetFromOptions(solver);
 
  auto dm = simple->getDM();
  auto X = createDMVector(dm);
  auto F = vectorDuplicate(X);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Solver start";
  CHKERR KSPSolve(solver, F, X);
  CHKERR VecGhostUpdateBegin(X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(X, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_REVERSE);
 
  MOFEM_LOG("INITIAL", Sev::inform)
      << "writing vector in binary to " << out_file_name << " ...";
  PetscViewer viewer;
  PetscViewerBinaryOpen(PETSC_COMM_WORLD, out_file_name, FILE_MODE_WRITE,
                        &viewer);
  VecView(X, viewer);
  PetscViewerDestroy(&viewer);
 
  MOFEM_LOG("INITIAL", Sev::inform) << "Solver done";
  MoFEMFunctionReturn(0);
}
//! [solveSystem]
 
MoFEMErrorCode PhotonDiffusion::outputResults() {
  MoFEMFunctionBegin;
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
 
  auto u_ptr = boost::make_shared<VectorDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE", u_ptr));
 
  using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(
 
          post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
          {{"PHOTON_FLUENCE_RATE", u_ptr}},
 
          {},
 
          {},
 
          {})
 
  );
 
  pipeline_mng->getDomainPostProcFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElementsPostProc();
  CHKERR post_proc_fe->writeFile("out_initial.h5m");
  MoFEMFunctionReturn(0);
}
 
//! [runProgram]
MoFEMErrorCode PhotonDiffusion::runProgram() {
  MoFEMFunctionBegin;
 
  CHKERR readMesh();
  CHKERR setupProblem();
  CHKERR initialCondition();
  CHKERR boundaryCondition();
  CHKERR assembleSystem();
  CHKERR solveSystem();
  if (output_volume)
    CHKERR outputResults();
  if (testing) {
    CHKERR createCommonData();
    CHKERR checkResults();
  }
  MoFEMFunctionReturn(0);
}
//! [runProgram]
 
//! [Check results]
MoFEMErrorCode PhotonDiffusion::checkResults() {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getOpDomainRhsPipeline().clear();
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateScalarFieldValues("PHOTON_FLUENCE_RATE",
                                       commonDataPtr->uAtPtsPtr));
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpError(commonDataPtr));
  CHKERR pipeline_mng->loopFiniteElements();
  CHKERR VecAssemblyBegin(commonDataPtr->L2Vec);
  CHKERR VecAssemblyEnd(commonDataPtr->L2Vec);
  CHKERR VecAssemblyBegin(commonDataPtr->resVec);
  CHKERR VecAssemblyEnd(commonDataPtr->resVec);
  double nrm2;
  CHKERR VecNorm(commonDataPtr->resVec, NORM_2, &nrm2);
  const double *array;
  CHKERR VecGetArrayRead(commonDataPtr->L2Vec, &array);
  if (mField.get_comm_rank() == 0)
    PetscPrintf(PETSC_COMM_SELF, "Error %6.4e Vec norm %6.4e\n",
                std::sqrt(array[0]), nrm2);
  CHKERR VecRestoreArrayRead(commonDataPtr->L2Vec, &array);
  constexpr double eps = 1e-8;
  if (nrm2 > eps)
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Not converged solution");
  MoFEMFunctionReturn(0);
}
//! [Check results]
 
int main(int argc, char *argv[]) {
 
  // Initialisation of MoFEM/PETSc and MOAB data structures
  const char param_file[] = "param_file.petsc";
  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(), "INITIAL"));
  LogManager::setLog("INITIAL");
  MOFEM_LOG_TAG("INITIAL", "initial_diffusion")
 
  // Error handling
  try {
    // Register MoFEM discrete manager in PETSc
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
 
    // Create MOAB instance
    moab::Core mb_instance;              // mesh database
    moab::Interface &moab = mb_instance; // mesh database interface
 
    // Create MoFEM instance
    MoFEM::Core core(moab);           // finite element database
    MoFEM::Interface &m_field = core; // finite element interface
 
    // Run the main analysis
    PhotonDiffusion heat_problem(m_field);
    CHKERR heat_problem.runProgram();
  }
  CATCH_ERRORS;
 
  // Finish work: cleaning memory, getting statistics, etc.
  MoFEM::Core::Finalize();
 
  return 0;
}