v0.14.0
lorentz_force.cpp

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

take in account higher order geometry

code panellisation, works serial only

Exaltation

\[ 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 \]

/**
* \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 <MoFEM.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, 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;
}
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
MoFEM::ForcesAndSourcesCore::getRuleHook
RuleHookFun getRuleHook
Hook to get rule.
Definition: ForcesAndSourcesCore.hpp:42
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::EntitiesFieldData::EntData::getFieldData
const VectorDouble & getFieldData() const
get dofs values
Definition: EntitiesFieldData.hpp:1254
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
FTensor::levi_civita
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
Definition: Levi_Civita.hpp:617
FTensor::Tensor2< double, 3, 3 >
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::invertTensor3by3
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:1588
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::FieldEvaluatorInterface::buildTree3D
MoFEMErrorCode buildTree3D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
Definition: FieldEvaluator.cpp:57
MoFEM::OpCalculateHcurlVectorCurl< 3, 3 >
Calculate curl of vector field.
Definition: UserDataOperators.hpp:2499
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::FieldEvaluatorInterface
Field evaluator interface.
Definition: FieldEvaluator.hpp:21
convert.type
type
Definition: convert.py:64
MoFEM::FieldEvaluatorInterface::evalFEAtThePoint3D
MoFEMErrorCode evalFEAtThePoint3D(const double *const point, const double distance, const std::string problem, const std::string finite_element, boost::shared_ptr< SetPtsData > data_ptr, int lower_rank, int upper_rank, boost::shared_ptr< CacheTuple > cache_ptr, MoFEMTypes bh=MF_EXIST, VERBOSITY_LEVELS verb=QUIET)
Evaluate field at artbitray position.
Definition: FieldEvaluator.cpp:364
main
int main(int argc, char *argv[])
Definition: lorentz_force.cpp:35
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:105
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::CoreInterface::cache_problem_entities
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
MoFEM::DMMoFEMCreateMoFEM
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: DMMoFEM.cpp:114
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
VolOp
VolEle::UserDataOperator VolOp
Definition: test_cache_on_entities.cpp:23
MoFEM::FieldEvaluatorInterface::getData
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
Definition: FieldEvaluator.hpp:107
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FTensor::Index< 'i', 3 >
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:22
convert.n
n
Definition: convert.py:82
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::DMMoFEMGetProblemPtr
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMoFEM.cpp:426
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
get_rule
auto get_rule
Definition: field_evaluator.cpp:73
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
dt
double dt
Definition: heat_method.cpp:26
MoFEM::Problem::getName
auto getName() const
Definition: ProblemsMultiIndices.hpp:372
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
MF_EXIST
@ MF_EXIST
Definition: definitions.h:113
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
help
static char help[]
Definition: lorentz_force.cpp:33