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.
\[ 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 \]
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
try {
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());
PetscBool flg_file;
CHKERR PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Lorenz force configure",
"none");
CHKERR PetscOptionsString(
"-my_file",
"mesh file name",
"",
"solution.h5m",
ierr = PetscOptionsEnd();
const char *option;
option = "";
DM dm;
CHKERR DMCreate(PETSC_COMM_WORLD, &dm);
CHKERR DMSetType(dm,
"DMMOFEM");
struct MyOpDebug :
public VolOp {
boost::shared_ptr<MatrixDouble> B;
MyOpDebug(decltype(B) b) :
VolOp(
"MAGNETIC_POTENTIAL", OPROW), B(b) {}
if (
type == MBEDGE && side == 0) {
std::cout << "found " << (*B) << endl;
}
}
};
const double dist = 1e-12;
const int nb_random_points = 5000;
const int nb_steps = 10000;
const int mod_step = 10;
const double velocity_scale = 0.1;
const double magnetic_field_scale = 1;
const double scale_box = 0.5;
auto vol_ele = data_at_pts->feMethodPtr.lock();
if (!vol_ele)
"Pointer to element does not exists");
auto get_rule = [&](
int order_row,
int order_col,
int order_data) {
return -1;
};
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,
auto &arrays_coord) {
ReadUtilIface *iface;
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](
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;
};
vector<double *> arrays_coord;
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 -
return t_rhs;
};
auto calc_lhs = [&](auto &t_B) {
for (auto ii : {0, 1, 2})
return t_lhs;
};
auto is_out = [&](auto &t_p) {
return true;
}
else if (t_p(
i) < bMin) {
return true;
}
return false;
};
auto cache_ptr = boost::make_shared<CacheTuple>();
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);
point.data(),
dist, prb_ptr->
getName(),
"MAGNETIC", data_at_pts,
if (B->size2())
for (int ii : {0, 1, 2})
t_B(ii) = (*B)(ii, 0);
else
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;
t_p(
i) = t_inv_lhs(
i,
j) * t_rhs(
j);
++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) {
CHKERR moab_charged_partices.write_file(
(
"step_" + boost::lexical_cast<std::string>(
t / mod_step) +
".vtk")
.c_str(),
"VTK", "");
}
}
std::cout << endl;
}
return 0;
}