v0.14.0
Loading...
Searching...
No Matches
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
31using namespace MoFEM;
32
33static char help[] = "...\n\n";
34
35int 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();
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();
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
86 using VolOp = VolumeElementForcesAndSourcesCore::UserDataOperator;
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]};
159 FTensor::Index<'i', 3> i;
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)};
180 FTensor::Index<'i', 3> i;
181 for (int n = 0; n != nb_random_points; ++n) {
182
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) {
214 FTensor::Index<'i', 3> i;
215 FTensor::Index<'j', 3> j;
216 FTensor::Index<'k', 3> k;
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) {
224 FTensor::Index<'i', 3> i;
225 FTensor::Index<'j', 3> j;
226 FTensor::Index<'k', 3> k;
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();
260 FTensor::Index<'i', 3> i;
261 FTensor::Index<'j', 3> j;
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 }
329
331
332 return 0;
333}
int main()
@ QUIET
#define CATCH_ERRORS
Catch errors.
@ MF_EXIST
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto get_rule
FTensor::Index< 'n', SPACE_DIM > n
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:483
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:118
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition DMMoFEM.cpp:412
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
FTensor::Index< 'i', SPACE_DIM > i
double dt
static char help[]
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
char mesh_file_name[255]
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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.
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
constexpr double t
plate stiffness
Definition plate.cpp:59
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:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorDouble & getFieldData() const
get dofs values
Field evaluator interface.
MoFEMErrorCode buildTree3D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
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.
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.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
RuleHookFun getRuleHook
Hook to get rule.
Calculate curl of vector field.
keeps basic data about problem
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolEle::UserDataOperator VolOp