MAX-1: Lorenz force

Note

Prerequisites of this tutorial include MAX-0: Magnetostatics

Note

Intended learning outcome:

• general structure of a program developed using MoFEM
• working with vectorial base fields like fields for H-curl space
• use of FieldEvaluatorInterface

/**
 * \file lorentz_force.cpp
 * \example lorentz_force.cpp
 * \brief Calculate Lorentz fore from magnetic field.
 *
 * It is not attemt to have accurate or realistic model of moving particles in
 * magnetic field. This example was create to show how field evaluator works
 * with precalculated magnetic field.
 *
 * \todo use physical quantities
 * \todo take in account higher order geometry
 * \todo code panellisation, works serial only
 *
 * Exaltation
 * \f[
 *  v_i = (p_{i+1} - p_{i-1}) / (2 \delta t) \\
 *  (p_{i-1} - 2 p_i + p_{i+1}) / \delta t^2 = \frac{q}{m} v_i \times B_i \\
 *  (p_{i-1} - 2 p_i + p_{i+1}) / \delta t^2 = \frac{q}{m} (p_{i+1} - p_{i-1}) /
 * (2 \delta t)  \times B_i \\
 *  (p_{i-1} - 2 p_i + p_{i+1}) / \delta t = \frac{q}{m} (p_{i+1} - p_{i-1})
 * \times B_i / 2  \\
 *  p_{i+1} / \delta t - p_{i+1} \times B_i / 2= (2 p_i - p_{i-1}) / \delta t -
 * p_{i-1} \times B_i / 2 \\ p_{i+1} (\mathbf{1} / \delta t - \frac{q}{m}
 * \mathbf{1} \times B_i / 2 )= (2 p_i - p_{i-1}) / \delta t - \frac{q}{m}
 * p_{i-1} \times B_i / 2 \f]
 *
 * \ingroup maxwell_element
 */
 
#include <BasicFiniteElements.hpp>
using namespace MoFEM;
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    moab::Core mb_instance;
    moab::Interface &moab = mb_instance;
 
    ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
    auto moab_comm_wrap =
        boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
    if (pcomm == NULL)
      pcomm =
          new ParallelComm(&moab, moab_comm_wrap->get_comm());
 
    // Read parameters from line command
    PetscBool flg_file;
    char mesh_file_name[255];
    CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Lorenz force configure",
                             "none");
    CHKERR PetscOptionsString("-my_file", "mesh file name", "", "solution.h5m",
                              mesh_file_name, 255, &flg_file);
    ierr = PetscOptionsEnd();
    CHKERRG(ierr);
 
    const char *option;
    option = "";
    CHKERR moab.load_file(mesh_file_name, 0, option);
 
    // Create mofem interface
    MoFEM::Core core(moab);
    MoFEM::Interface &m_field = core;
 
    CHKERR m_field.build_fields();
    CHKERR m_field.build_finite_elements();
    CHKERR m_field.build_adjacencies(BitRefLevel().set(0));
 
    // set up DM
    DM dm;
    CHKERR DMRegister_MoFEM("DMMOFEM");
    CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
    CHKERR DMSetType(dm, "DMMOFEM");
    CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "MAGNETIC_PROBLEM",
                              BitRefLevel().set(0));
    CHKERR DMSetFromOptions(dm);
    // add elements to blockData.dM
    CHKERR DMMoFEMAddElement(dm, "MAGNETIC");
    CHKERR DMSetUp(dm);
 
    using VolEle = VolumeElementForcesAndSourcesCore;
    using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
    // using SetPtsData = FieldEvaluatorInterface::SetPtsData;
    // using SetPts = FieldEvaluatorInterface::SetPts;
 
    /**
     * @brief Only for debuging
     */
    struct MyOpDebug : public VolOp {
 
      boost::shared_ptr<MatrixDouble> B;
      MyOpDebug(decltype(B) b) : VolOp("MAGNETIC_POTENTIAL", OPROW), B(b) {}
 
      MoFEMErrorCode doWork(int side, EntityType type,
                            EntitiesFieldData::EntData &data) {
        MoFEMFunctionBegin;
        if (type == MBEDGE && side == 0) {
          std::cout << "found " << (*B) << endl;
          std::cout << "data " << data.getFieldData() << endl;
        }
 
        MoFEMFunctionReturn(0);
      }
    };
 
    const double dist = 1e-12;             // Distance for tree search
    const int nb_random_points = 5000;     // number of points
    const int nb_steps = 10000;            // number of time steps
    const int mod_step = 10;               // save every step
    const double dt = 1e-5;                // Time step size
    const double velocity_scale = 0.1;     // scale velocity vector
    const double magnetic_field_scale = 1; // scale magnetic field vector
    const double scale_box = 0.5; // scale box where partices are placed
 
    FieldEvaluatorInterface *field_eval_ptr;
    CHKERR m_field.getInterface(field_eval_ptr);
 
    // Get access to data
    auto data_at_pts = field_eval_ptr->getData<VolEle>();
    auto vol_ele = data_at_pts->feMethodPtr.lock();
    if (!vol_ele)
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "Pointer to element does not exists");
    auto get_rule = [&](int order_row, int order_col, int order_data) {
      return -1;
    };
    vol_ele->getRuleHook = get_rule;
 
    boost::shared_ptr<MatrixDouble> B(new MatrixDouble());
    vol_ele->getOpPtrVector().push_back(
        new OpCalculateHcurlVectorCurl<3, 3>("MAGNETIC_POTENTIAL", B));
    // vol_ele->getOpPtrVector().push_back(new MyOpDebug(B));
 
    const MoFEM::Problem *prb_ptr;
    CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
 
    CHKERR field_eval_ptr->buildTree3D(data_at_pts, "MAGNETIC");
    BoundBox box;
    CHKERR data_at_pts->treePtr->get_bounding_box(box);
 
    const double bMin = box.bMin[0];
    const double bMax = box.bMax[0];
 
    auto create_vertices = [nb_random_points, bMin, bMax,
                            scale_box](moab::Interface &moab, auto &verts,
                                       auto &arrays_coord) {
      MoFEMFunctionBegin;
      ReadUtilIface *iface;
      EntityHandle startv;
      CHKERR moab.query_interface(iface);
      CHKERR iface->get_node_coords(3, nb_random_points, 0, startv,
                                    arrays_coord);
      FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords = {
          arrays_coord[0], arrays_coord[1], arrays_coord[2]};
      FTensor::Index<'i', 3> i;
      verts = Range(startv, startv + nb_random_points - 1);
      for (int n = 0; n != nb_random_points; ++n) {
        t_coords(0) = 0;
        for (auto ii : {1, 2}) {
          t_coords(ii) = scale_box * (bMax - bMin) *
                             (std::rand() / static_cast<double>(RAND_MAX)) -
                         bMax * scale_box;
        }
        ++t_coords;
      }
      MoFEMFunctionReturn(0);
    };
 
    auto set_positions = [nb_random_points, dt, velocity_scale](
                             moab::Interface &moab, auto &arrays_coord) {
      MatrixDouble init_pos(3, nb_random_points);
      FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_coords = {
          arrays_coord[0], arrays_coord[1], arrays_coord[2]};
      FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_init_coords = {
          &init_pos(0, 0), &init_pos(1, 0), &init_pos(2, 0)};
      FTensor::Index<'i', 3> i;
      for (int n = 0; n != nb_random_points; ++n) {
 
        FTensor::Tensor1<double, 3> t_velocity;
        for (auto ii : {0, 1, 2})
          t_velocity(ii) = (rand() / static_cast<double>(RAND_MAX) - 0.5);
        t_velocity(i) /= std::sqrt(t_velocity(i) * t_velocity(i));
        t_init_coords(i) = t_coords(i) + dt * velocity_scale * t_velocity(i);
 
        ++t_coords;
        ++t_init_coords;
      }
      return init_pos;
    };
 
    moab::Core mb_charged_partices;
    moab::Interface &moab_charged_partices = mb_charged_partices;
    vector<double *> arrays_coord;
    Range verts;
    CHKERR create_vertices(moab_charged_partices, verts, arrays_coord);
    auto init_positions = set_positions(moab, arrays_coord);
 
    auto get_t_coords = [&]() {
      return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
          arrays_coord[0], arrays_coord[1], arrays_coord[2]);
    };
 
    auto get_t_init_coords = [&]() {
      return FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3>(
          &init_positions(0, 0), &init_positions(1, 0), &init_positions(2, 0));
    };
 
    auto calc_rhs = [&](auto &t_p, auto &t_init_p, auto &t_B) {
      FTensor::Tensor1<double, 3> t_rhs;
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'k', 3> k;
      t_rhs(k) = (2 * t_p(k) - t_init_p(k)) / dt -
                 levi_civita(j, i, k) * t_init_p(i) * t_B(j);
      return t_rhs;
    };
 
    auto calc_lhs = [&](auto &t_B) {
      FTensor::Tensor2<double, 3, 3> t_lhs;
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
      FTensor::Index<'k', 3> k;
      t_lhs(i, k) = levi_civita(j, i, k) * (-t_B(j));
      for (auto ii : {0, 1, 2})
        t_lhs(ii, ii) += 1 / dt;
      return t_lhs;
    };
 
    // auto set_periodicity = [&](auto &t_p, auto &t_init_p) {
    //   for (int i : {0, 1, 2})
    //     if (t_p(i) > bMax) {
    //       t_p(i) -= 2 * bMax;
    //       t_init_p(i) -= 2 * bMax;
    //     } else if (t_p(i) < bMin) {
    //       t_p(i) -= 2 * bMin;
    //       t_init_p(i) -= 2 * bMin;
    //     }
    // };
 
    auto is_out = [&](auto &t_p) {
      for (int i : {0, 1, 2})
        if (t_p(i) > bMax) {
          return true;
        } else if (t_p(i) < bMin) {
          return true;
        }
      return false;
    };
 
    auto cache_ptr = boost::make_shared<CacheTuple>();
    CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
 
    auto calc_position = [&]() {
      auto t_p = get_t_coords();
      auto t_init_p = get_t_init_coords();
      FTensor::Index<'i', 3> i;
      FTensor::Index<'j', 3> j;
 
      for (int n = 0; n != nb_random_points; ++n) {
 
        if (is_out(t_p)) {
          ++t_p;
          ++t_init_p;
          continue;
        }
 
        std::array<double, 3> point = {t_p(0), t_p(1), t_p(2)};
        data_at_pts->setEvalPoints(point.data(), 1);
 
        CHKERR field_eval_ptr->evalFEAtThePoint3D(
            point.data(), dist, prb_ptr->getName(), "MAGNETIC", data_at_pts,
            m_field.get_comm_rank(), m_field.get_comm_rank(), cache_ptr,
            MF_EXIST, QUIET);
 
        FTensor::Tensor1<double, 3> t_B;
 
        if (B->size2())
          for (int ii : {0, 1, 2})
            t_B(ii) = (*B)(ii, 0);
        else
          t_B(i) = 0;
 
        t_B(i) *= magnetic_field_scale * 0.5;
 
        auto t_rhs = calc_rhs(t_p, t_init_p, t_B);
        auto t_lhs = calc_lhs(t_B);
 
        double det;
        CHKERR determinantTensor3by3(t_lhs, det);
        FTensor::Tensor2<double, 3, 3> t_inv_lhs;
        CHKERR invertTensor3by3(t_lhs, det, t_inv_lhs);
 
        t_init_p(i) = t_p(i);
        t_p(i) = t_inv_lhs(i, j) * t_rhs(j);
 
        // set_periodicity(t_p, t_init_p);
 
        ++t_p;
        ++t_init_p;
      }
    };
 
    for (int t = 0; t != nb_steps; ++t) {
 
      std::string print_step =
          "Step : " + boost::lexical_cast<std::string>(t) + "\r";
      std::cout << print_step << std::flush;
 
      calc_position();
 
      if ((t % mod_step) == 0) {
        // write points
        CHKERR moab_charged_partices.write_file(
            ("step_" + boost::lexical_cast<std::string>(t / mod_step) + ".vtk")
                .c_str(),
            "VTK", "");
      }
    }
 
    std::cout << endl;
 
    CHKERR DMDestroy(&dm);
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return 0;
}