v0.13.0
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
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#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();
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();
// set up DM
DM dm;
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 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,
if (type == MBEDGE && side == 0) {
std::cout << "found " << (*B) << endl;
std::cout << "data " << data.getFieldData() << endl;
}
}
};
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>("MAGNETIC_POTENTIAL", B));
// vol_ele->getOpPtrVector().push_back(new MyOpDebug(B));
const MoFEM::Problem *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) {
ReadUtilIface *iface;
EntityHandle startv;
CHKERR moab.query_interface(iface);
CHKERR iface->get_node_coords(3, nb_random_points, 0, startv,
arrays_coord);
arrays_coord[0], arrays_coord[1], arrays_coord[2]};
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;
}
};
auto set_positions = [nb_random_points, dt, velocity_scale](
moab::Interface &moab, auto &arrays_coord) {
MatrixDouble init_pos(3, nb_random_points);
arrays_coord[0], arrays_coord[1], arrays_coord[2]};
&init_pos(0, 0), &init_pos(1, 0), &init_pos(2, 0)};
for (int n = 0; n != nb_random_points; ++n) {
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 = [&]() {
arrays_coord[0], arrays_coord[1], arrays_coord[2]);
};
auto get_t_init_coords = [&]() {
&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) {
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) {
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();
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,
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 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);
}
return 0;
}
static Index< 'n', 3 > n
int main(int argc, char *argv[])
static char help[]
EntitiesFieldData::EntData EntData
@ QUIET
Definition: definitions.h:221
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_EXIST
Definition: definitions.h:113
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto get_rule
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:130
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:461
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:384
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
FTensor::Index< 'i', SPACE_DIM > i
double dt
Definition: heat_method.cpp:40
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
char mesh_file_name[255]
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1172
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1161
constexpr double t
plate stiffness
Definition: plate.cpp:76
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
RuleHookFun getRuleHook
Hook to get rule.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator