v0.14.0
lorentz_force.cpp
Go to the documentation of this file.
1 /**
2  * \file lorentz_force.cpp
3  * \example lorentz_force.cpp
4  * \brief Calculate Lorentz fore from magnetic field.
5  *
6  * It is not attemt to have accurate or realistic model of moving particles in
7  * magnetic field. This example was create to show how field evaluator works
8  * with precalculated magnetic field.
9  *
10  * \todo use physical quantities
11  * \todo take in account higher order geometry
12  * \todo code panellisation, works serial only
13  *
14  * Exaltation
15  * \f[
16  * v_i = (p_{i+1} - p_{i-1}) / (2 \delta t) \\
17  * (p_{i-1} - 2 p_i + p_{i+1}) / \delta t^2 = \frac{q}{m} v_i \times B_i \\
18  * (p_{i-1} - 2 p_i + p_{i+1}) / \delta t^2 = \frac{q}{m} (p_{i+1} - p_{i-1}) /
19  * (2 \delta t) \times B_i \\
20  * (p_{i-1} - 2 p_i + p_{i+1}) / \delta t = \frac{q}{m} (p_{i+1} - p_{i-1})
21  * \times B_i / 2 \\
22  * p_{i+1} / \delta t - p_{i+1} \times B_i / 2= (2 p_i - p_{i-1}) / \delta t -
23  * p_{i-1} \times B_i / 2 \\ p_{i+1} (\mathbf{1} / \delta t - \frac{q}{m}
24  * \mathbf{1} \times B_i / 2 )= (2 p_i - p_{i-1}) / \delta t - \frac{q}{m}
25  * p_{i-1} \times B_i / 2 \f]
26  *
27  * \ingroup maxwell_element
28  */
29 
30 #include <BasicFiniteElements.hpp>
31 using namespace MoFEM;
32 
33 static char help[] = "...\n\n";
34 
35 int main(int argc, char *argv[]) {
36 
37  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
38 
39  try {
40 
41  moab::Core mb_instance;
42  moab::Interface &moab = mb_instance;
43 
44  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
45  auto moab_comm_wrap =
46  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
47  if (pcomm == NULL)
48  pcomm =
49  new ParallelComm(&moab, moab_comm_wrap->get_comm());
50 
51  // Read parameters from line command
52  PetscBool flg_file;
53  char mesh_file_name[255];
54  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Lorenz force configure",
55  "none");
56  CHKERR PetscOptionsString("-my_file", "mesh file name", "", "solution.h5m",
57  mesh_file_name, 255, &flg_file);
58  ierr = PetscOptionsEnd();
59  CHKERRG(ierr);
60 
61  const char *option;
62  option = "";
63  CHKERR moab.load_file(mesh_file_name, 0, option);
64 
65  // Create mofem interface
66  MoFEM::Core core(moab);
67  MoFEM::Interface &m_field = core;
68 
69  CHKERR m_field.build_fields();
70  CHKERR m_field.build_finite_elements();
71  CHKERR m_field.build_adjacencies(BitRefLevel().set(0));
72 
73  // set up DM
74  DM dm;
75  CHKERR DMRegister_MoFEM("DMMOFEM");
76  CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
77  CHKERR DMSetType(dm, "DMMOFEM");
78  CHKERR DMMoFEMCreateMoFEM(dm, &m_field, "MAGNETIC_PROBLEM",
79  BitRefLevel().set(0));
80  CHKERR DMSetFromOptions(dm);
81  // add elements to blockData.dM
82  CHKERR DMMoFEMAddElement(dm, "MAGNETIC");
83  CHKERR DMSetUp(dm);
84 
87  // using SetPtsData = FieldEvaluatorInterface::SetPtsData;
88  // using SetPts = FieldEvaluatorInterface::SetPts;
89 
90  /**
91  * @brief Only for debuging
92  */
93  struct MyOpDebug : public VolOp {
94 
95  boost::shared_ptr<MatrixDouble> B;
96  MyOpDebug(decltype(B) b) : VolOp("MAGNETIC_POTENTIAL", OPROW), B(b) {}
97 
98  MoFEMErrorCode doWork(int side, EntityType type,
101  if (type == MBEDGE && side == 0) {
102  std::cout << "found " << (*B) << endl;
103  std::cout << "data " << data.getFieldData() << endl;
104  }
105 
107  }
108  };
109 
110  const double dist = 1e-12; // Distance for tree search
111  const int nb_random_points = 5000; // number of points
112  const int nb_steps = 10000; // number of time steps
113  const int mod_step = 10; // save every step
114  const double dt = 1e-5; // Time step size
115  const double velocity_scale = 0.1; // scale velocity vector
116  const double magnetic_field_scale = 1; // scale magnetic field vector
117  const double scale_box = 0.5; // scale box where partices are placed
118 
119  FieldEvaluatorInterface *field_eval_ptr;
120  CHKERR m_field.getInterface(field_eval_ptr);
121 
122  // Get access to data
123  auto data_at_pts = field_eval_ptr->getData<VolEle>();
124  auto vol_ele = data_at_pts->feMethodPtr.lock();
125  if (!vol_ele)
126  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
127  "Pointer to element does not exists");
128  auto get_rule = [&](int order_row, int order_col, int order_data) {
129  return -1;
130  };
131  vol_ele->getRuleHook = get_rule;
132 
133  boost::shared_ptr<MatrixDouble> B(new MatrixDouble());
134  vol_ele->getOpPtrVector().push_back(
135  new OpCalculateHcurlVectorCurl<3, 3>("MAGNETIC_POTENTIAL", B));
136  // vol_ele->getOpPtrVector().push_back(new MyOpDebug(B));
137 
138  const MoFEM::Problem *prb_ptr;
139  CHKERR DMMoFEMGetProblemPtr(dm, &prb_ptr);
140 
141  CHKERR field_eval_ptr->buildTree3D(data_at_pts, "MAGNETIC");
142  BoundBox box;
143  CHKERR data_at_pts->treePtr->get_bounding_box(box);
144 
145  const double bMin = box.bMin[0];
146  const double bMax = box.bMax[0];
147 
148  auto create_vertices = [nb_random_points, bMin, bMax,
149  scale_box](moab::Interface &moab, auto &verts,
150  auto &arrays_coord) {
152  ReadUtilIface *iface;
153  EntityHandle startv;
154  CHKERR moab.query_interface(iface);
155  CHKERR iface->get_node_coords(3, nb_random_points, 0, startv,
156  arrays_coord);
158  arrays_coord[0], arrays_coord[1], arrays_coord[2]};
160  verts = Range(startv, startv + nb_random_points - 1);
161  for (int n = 0; n != nb_random_points; ++n) {
162  t_coords(0) = 0;
163  for (auto ii : {1, 2}) {
164  t_coords(ii) = scale_box * (bMax - bMin) *
165  (std::rand() / static_cast<double>(RAND_MAX)) -
166  bMax * scale_box;
167  }
168  ++t_coords;
169  }
171  };
172 
173  auto set_positions = [nb_random_points, dt, velocity_scale](
174  moab::Interface &moab, auto &arrays_coord) {
175  MatrixDouble init_pos(3, nb_random_points);
177  arrays_coord[0], arrays_coord[1], arrays_coord[2]};
179  &init_pos(0, 0), &init_pos(1, 0), &init_pos(2, 0)};
181  for (int n = 0; n != nb_random_points; ++n) {
182 
183  FTensor::Tensor1<double, 3> t_velocity;
184  for (auto ii : {0, 1, 2})
185  t_velocity(ii) = (rand() / static_cast<double>(RAND_MAX) - 0.5);
186  t_velocity(i) /= std::sqrt(t_velocity(i) * t_velocity(i));
187  t_init_coords(i) = t_coords(i) + dt * velocity_scale * t_velocity(i);
188 
189  ++t_coords;
190  ++t_init_coords;
191  }
192  return init_pos;
193  };
194 
195  moab::Core mb_charged_partices;
196  moab::Interface &moab_charged_partices = mb_charged_partices;
197  vector<double *> arrays_coord;
198  Range verts;
199  CHKERR create_vertices(moab_charged_partices, verts, arrays_coord);
200  auto init_positions = set_positions(moab, arrays_coord);
201 
202  auto get_t_coords = [&]() {
204  arrays_coord[0], arrays_coord[1], arrays_coord[2]);
205  };
206 
207  auto get_t_init_coords = [&]() {
209  &init_positions(0, 0), &init_positions(1, 0), &init_positions(2, 0));
210  };
211 
212  auto calc_rhs = [&](auto &t_p, auto &t_init_p, auto &t_B) {
217  t_rhs(k) = (2 * t_p(k) - t_init_p(k)) / dt -
218  levi_civita(j, i, k) * t_init_p(i) * t_B(j);
219  return t_rhs;
220  };
221 
222  auto calc_lhs = [&](auto &t_B) {
227  t_lhs(i, k) = levi_civita(j, i, k) * (-t_B(j));
228  for (auto ii : {0, 1, 2})
229  t_lhs(ii, ii) += 1 / dt;
230  return t_lhs;
231  };
232 
233  // auto set_periodicity = [&](auto &t_p, auto &t_init_p) {
234  // for (int i : {0, 1, 2})
235  // if (t_p(i) > bMax) {
236  // t_p(i) -= 2 * bMax;
237  // t_init_p(i) -= 2 * bMax;
238  // } else if (t_p(i) < bMin) {
239  // t_p(i) -= 2 * bMin;
240  // t_init_p(i) -= 2 * bMin;
241  // }
242  // };
243 
244  auto is_out = [&](auto &t_p) {
245  for (int i : {0, 1, 2})
246  if (t_p(i) > bMax) {
247  return true;
248  } else if (t_p(i) < bMin) {
249  return true;
250  }
251  return false;
252  };
253 
254  auto cache_ptr = boost::make_shared<CacheTuple>();
255  CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
256 
257  auto calc_position = [&]() {
258  auto t_p = get_t_coords();
259  auto t_init_p = get_t_init_coords();
262 
263  for (int n = 0; n != nb_random_points; ++n) {
264 
265  if (is_out(t_p)) {
266  ++t_p;
267  ++t_init_p;
268  continue;
269  }
270 
271  std::array<double, 3> point = {t_p(0), t_p(1), t_p(2)};
272  data_at_pts->setEvalPoints(point.data(), 1);
273 
274  CHKERR field_eval_ptr->evalFEAtThePoint3D(
275  point.data(), dist, prb_ptr->getName(), "MAGNETIC", data_at_pts,
276  m_field.get_comm_rank(), m_field.get_comm_rank(), cache_ptr,
277  MF_EXIST, QUIET);
278 
280 
281  if (B->size2())
282  for (int ii : {0, 1, 2})
283  t_B(ii) = (*B)(ii, 0);
284  else
285  t_B(i) = 0;
286 
287  t_B(i) *= magnetic_field_scale * 0.5;
288 
289  auto t_rhs = calc_rhs(t_p, t_init_p, t_B);
290  auto t_lhs = calc_lhs(t_B);
291 
292  double det;
293  CHKERR determinantTensor3by3(t_lhs, det);
295  CHKERR invertTensor3by3(t_lhs, det, t_inv_lhs);
296 
297  t_init_p(i) = t_p(i);
298  t_p(i) = t_inv_lhs(i, j) * t_rhs(j);
299 
300  // set_periodicity(t_p, t_init_p);
301 
302  ++t_p;
303  ++t_init_p;
304  }
305  };
306 
307  for (int t = 0; t != nb_steps; ++t) {
308 
309  std::string print_step =
310  "Step : " + boost::lexical_cast<std::string>(t) + "\r";
311  std::cout << print_step << std::flush;
312 
313  calc_position();
314 
315  if ((t % mod_step) == 0) {
316  // write points
317  CHKERR moab_charged_partices.write_file(
318  ("step_" + boost::lexical_cast<std::string>(t / mod_step) + ".vtk")
319  .c_str(),
320  "VTK", "");
321  }
322  }
323 
324  std::cout << endl;
325 
326  CHKERR DMDestroy(&dm);
327  }
328  CATCH_ERRORS;
329 
331 
332  return 0;
333 }
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:1244
BasicFiniteElements.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:1566
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
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:59
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:23
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:1518
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:1140
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