WarpX
Loading...
Searching...
No Matches
MatrixPC.H
Go to the documentation of this file.
1/* Copyright 2024 Debojyoti Ghosh
2 *
3 * This file is part of WarpX.
4 *
5 * License: BSD-3-Clause-LBNL
6 */
7#ifndef MATRIX_PC_H_
8#define MATRIX_PC_H_
9
10#include "Fields.H"
11#include "Utils/WarpXConst.H"
12#include "Utils/TextMsg.H"
13#include "Preconditioner.H"
14
15#include <ablastr/fields/MultiFabRegister.H>
16
17#include <AMReX.H>
18#include <AMReX_ParmParse.H>
19#include <AMReX_Array.H>
20#include <AMReX_Vector.H>
21#include <AMReX_MultiFab.H>
22
23namespace MatrixPCUtils
24{
25 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
26 bool insertOrAdd( const int a_cidx,
27 const amrex::Real a_val,
28 int* const a_cidxs,
29 amrex::Real* const a_aij,
30 const int a_nnz,
31 int& a_ncol )
32 {
33 if (a_cidx < 0) { return true; /* outside domain */ }
34 int loc = -1;
35 for (int icol = 0; icol < std::min(a_ncol,a_nnz); icol++) {
36 if (a_cidxs[icol] == a_cidx) {
37 loc = icol;
38 break;
39 }
40 }
41 if (loc < 0) {
42 a_ncol++;
43 if (a_ncol > a_nnz) { return false; }
44 else {
45 // column index not found; add new entry
46 a_cidxs[a_ncol-1] = a_cidx;
47 a_aij[a_ncol-1] = a_val;
48 }
49 } else {
50 // column index already exists; add to value
51 a_aij[loc] += a_val;
52 }
53 return true;
54 }
55}
56
74
75template <class T, class Ops>
76class MatrixPC : public Preconditioner<T,Ops>
77{
78 public:
79
80 using RT = typename T::value_type;
81
85 MatrixPC () = default;
86
90 ~MatrixPC () override = default;
91
92 // Prohibit move and copy operations
93 MatrixPC(const MatrixPC&) = delete;
94 MatrixPC& operator=(const MatrixPC&) = delete;
95 MatrixPC(MatrixPC&&) noexcept = delete;
96 MatrixPC& operator=(MatrixPC&&) noexcept = delete;
97
101 void Define (const T&, Ops* const) override;
102
106 void Update (const T& a_U) override;
107
116 int Assemble (const T& a_U);
117
126 void Apply (T&, const T&) override;
127
128 inline void getPCMatrix (amrex::Gpu::DeviceVector<int>& a_r_indices_g,
129 amrex::Gpu::DeviceVector<int>& a_num_nz,
130 amrex::Gpu::DeviceVector<int>& a_c_indices_g,
131 amrex::Gpu::DeviceVector<RT>& a_a_ij,
132 int& a_n, int& a_ncols_max ) override
133 {
134 a_n = m_ndofs_l;
135 a_ncols_max = m_pc_mat_nnz;
136
137 a_r_indices_g.resize(m_r_indices_g.size());
138 amrex::Gpu::copyAsync( amrex::Gpu::deviceToDevice,
139 m_r_indices_g.begin(),
140 m_r_indices_g.end(),
141 a_r_indices_g.begin() );
142
143 a_num_nz.resize(m_num_nz.size());
144 amrex::Gpu::copyAsync( amrex::Gpu::deviceToDevice,
145 m_num_nz.begin(),
146 m_num_nz.end(),
147 a_num_nz.begin() );
148
149 a_c_indices_g.resize(m_c_indices_g.size());
150 amrex::Gpu::copyAsync( amrex::Gpu::deviceToDevice,
151 m_c_indices_g.begin(),
152 m_c_indices_g.end(),
153 a_c_indices_g.begin() );
154
155 a_a_ij.resize(m_a_ij.size());
156 amrex::Gpu::copyAsync( amrex::Gpu::deviceToDevice,
157 m_a_ij.begin(),
158 m_a_ij.end(),
159 a_a_ij.begin() );
160 amrex::Gpu::streamSynchronize();
161 }
162
166 void printParameters() const override;
167
171 [[nodiscard]] inline bool IsDefined () const override { return m_is_defined; }
172
173 inline void setName (const std::string& a_name) override { m_name = a_name; }
174
175 protected:
176
177 bool m_is_defined = false;
178 bool m_verbose = true;
179
180 Ops* m_ops = nullptr;
181
182 int m_num_amr_levels = 0;
183 amrex::Vector<amrex::Geometry> m_geom;
184
185 int m_ndofs_l = 0;
186 int m_ndofs_g = 0;
187 bool m_pc_diag_only = false;
188 int m_pc_mat_nnz = 1;
189 bool m_include_mass_matrices = false;
190
191 std::string m_name = "noname";
192
193 amrex::Gpu::DeviceVector<int> m_r_indices_g;
194 amrex::Gpu::DeviceVector<int> m_num_nz;
195 amrex::Gpu::DeviceVector<int> m_c_indices_g;
196 amrex::Gpu::DeviceVector<amrex::Real> m_a_ij;
197
198 const amrex::Vector<amrex::Array<amrex::MultiFab*,3>>* m_bcoefs = nullptr;
199
200 int m_num_realloc = 0;
201
205 void readParameters();
206
207 private:
208
209};
210
211template <class T, class Ops>
212void MatrixPC<T,Ops>::printParameters() const
213{
214 using namespace amrex;
215 Print() << m_name << " verbose: " << (m_verbose?"true":"false") << "\n";
216 Print() << m_name << " pc_diagonal_only: " << (m_pc_diag_only?"true":"false") << "\n";
217 Print() << m_name << " include_mass_matrices: " << (m_include_mass_matrices?"true":"false") << "\n";
218}
219
220template <class T, class Ops>
221void MatrixPC<T,Ops>::readParameters()
222{
223 const amrex::ParmParse pp(m_name);
224 pp.query("verbose", m_verbose);
225 pp.query("pc_diagonal_only", m_pc_diag_only);
226}
227
228template <class T, class Ops>
229void MatrixPC<T,Ops>::Define ( const T& a_U,
230 Ops* const a_ops )
231{
232 BL_PROFILE("MatrixPC::Define()");
233 using namespace amrex;
234
235 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
236 !IsDefined(),
237 "MatrixPC::Define() called on defined object" );
238 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
239 (a_ops != nullptr),
240 "MatrixPC::Define(): a_ops is nullptr" );
241 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
242 a_U.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
243 "MatrixPC::Define() must be called with Efield_fp type");
244
245 m_ops = a_ops;
246 // read preconditioner parameters
247 readParameters();
248
249 // a_U is not needed
250 amrex::ignore_unused(a_U);
251
252 // Set number of AMR levels and create geometry, grids, and
253 // distribution mapping vectors.
254 m_num_amr_levels = m_ops->numAMRLevels();
255 if (m_num_amr_levels > 1) {
256 WARPX_ABORT_WITH_MESSAGE("MatrixPC::Define(): m_num_amr_levels > 1");
257 }
258 m_geom.resize(m_num_amr_levels);
259 for (int n = 0; n < m_num_amr_levels; n++) {
260 m_geom[n] = m_ops->GetGeometry(n);
261 }
262
263 m_ndofs_l = a_U.nDOF_local();
264 m_ndofs_g = a_U.nDOF_global();
265
266 auto n_rows = size_t(m_ndofs_l);
267 auto n_cols = size_t(m_pc_mat_nnz) * size_t(m_ndofs_l);
268
269 m_r_indices_g.resize(n_rows);
270 m_num_nz.resize(n_rows);
271 m_c_indices_g.resize(n_cols);
272 m_a_ij.resize(n_cols);
273
274
275 m_bcoefs = m_ops->GetMassMatricesCoeff();
276 if (m_bcoefs != nullptr) {
277 m_include_mass_matrices = true;
278 }
279
280 m_is_defined = true;
281}
282
283template <class T, class Ops>
284void MatrixPC<T,Ops>::Update (const T& a_U)
285{
286 BL_PROFILE("MatrixPC::Update()");
287 using namespace amrex;
288
289 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
290 IsDefined(),
291 "MatrixPC::Update() called on undefined object" );
292
293 while(1) {
294
295 auto nnz_diff = Assemble(a_U);
296 AMREX_ALWAYS_ASSERT(nnz_diff >= 0);
297 if (nnz_diff) {
298
299 m_pc_mat_nnz += nnz_diff;
300 m_num_realloc++;
301
302 } else {
303 break;
304 }
305 }
306
307 if (m_num_realloc > 1) {
308 std::stringstream warning_message;
309 warning_message << "Number of times arrays were reallocated due to new nonzero elements "
310 << "is greater than 1 (" << m_num_realloc <<"). This is unexpected.\n";
311 ablastr::warn_manager::WMRecordWarning("MatrixPC", warning_message.str());
312 }
313
314}
315
316template <class T, class Ops>
317int MatrixPC<T,Ops>::Assemble (const T& a_U)
318{
319 // Assemble the sparse matrix representation of the preconditioner
320 // A = curl (alpha * curl []) + M
321 // where M is the mass matrix. The following data is set in this function:
322 // - m_r_indices_g: integer array of size n with the global row indices
323 // - m_num_nz: integer array of size n with the number of non-zero elements
324 // in each row
325 // - m_c_indices: integer array of size n*ncmax with the global column indices
326 // of non-zero elements in each row (row-major)
327 // - m_a_ij: real-type array of size n*ncmax with the matrix element values
328 // (row-major format)
329 // where n is the local number of rows, and ncmax is the maximum number of non-zero
330 // elements per row.
331
332 BL_PROFILE("MatrixPC::Assemble()");
333 using namespace amrex;
334
335 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
336 IsDefined(),
337 "MatrixPC::Assemble() called on undefined object" );
338
339 // set the alpha coefficient for the curl-curl op
340 const RT thetaDt = m_ops->GetThetaForPC()*this->m_dt;
341 const RT alpha = (thetaDt*PhysConst::c) * (thetaDt*PhysConst::c);
342 if (m_verbose) {
343 amrex::Print() << "Updating " << m_name
344 << ": theta*dt = " << thetaDt << ", "
345 << " coefficients: "
346 << "alpha = " << alpha << "\n";
347 }
348
349 // Get DOF object from a_U
350 const auto& dofs_obj = a_U.getDOFsObject();
351 const auto& dofs_mfarrvec = dofs_obj->m_array;
352 AMREX_ALWAYS_ASSERT(m_ndofs_l == dofs_obj->m_nDoFs_l);
353 AMREX_ALWAYS_ASSERT(m_ndofs_g == dofs_obj->m_nDoFs_g);
354
355 m_r_indices_g.clear();
356 m_num_nz.clear();
357 m_c_indices_g.clear();
358 m_a_ij.clear();
359
360 auto n_rows = size_t(m_ndofs_l);
361 auto n_cols = size_t(m_pc_mat_nnz) * size_t(m_ndofs_l);
362
363 m_r_indices_g.resize(n_rows);
364 m_num_nz.resize(n_rows);
365 m_c_indices_g.resize(n_cols);
366 m_a_ij.resize(n_cols);
367
368 auto r_indices_g_ptr = m_r_indices_g.data();
369 auto num_nz_ptr = m_num_nz.data();
370 auto c_indices_g_ptr = m_c_indices_g.data();
371 auto a_ij_ptr = m_a_ij.data();
372
373 const auto nnz_max = m_pc_mat_nnz;
374 auto nnz_actual = nnz_max;
375
376 for (int lev = 0; lev < m_num_amr_levels; lev++) {
377
378 auto ncomp = dofs_mfarrvec[lev][0]->nComp();
379 AMREX_ALWAYS_ASSERT(ncomp == 2); // local, global
380
381 const auto& geom = m_geom[lev];
382 const auto dxi = geom.InvCellSizeArray();
383
384 Gpu::Buffer<int> nnz_actual_d({nnz_max});
385 auto* nnz_actual_ptr = nnz_actual_d.data();
386
387 for (int dir = 0; dir < 3; dir++) {
388
389 const amrex::MultiFab* BC_mask_Edir = m_ops->GetCurl2BCmask(lev,dir);
390
391 for (amrex::MFIter mfi(*dofs_mfarrvec[lev][dir]); mfi.isValid(); ++mfi) {
392
393 const amrex::Box bx = mfi.tilebox();
394 const amrex::Box full_bx = mfi.fabbox();
395
396 // BC mask values for in-plane components of curl curl E operator for 2D
397 const auto BC_mask_Edir_arr = BC_mask_Edir->const_array(mfi);
398
399 auto dof_arr = dofs_mfarrvec[lev][dir]->const_array(mfi);
400
401#if defined(WARPX_DIM_RCYLINDER)
402 amrex::ignore_unused(dxi, BC_mask_Edir_arr);
403 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RCYLINDER");
404#elif defined(WARPX_DIM_RSPHERE)
405 amrex::ignore_unused(dxi, BC_mask_Edir_arr);
406 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RSPHERE");
407#elif defined(WARPX_DIM_RZ)
408 amrex::ignore_unused(dxi, BC_mask_Edir_arr);
409 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RZ");
410#elif defined(WARPX_DIM_XZ)
411 int tdir = -1;
412 if (dir == 0) { tdir = 2; }
413 else if (dir == 2) { tdir = 0; }
414 else { tdir = 1; }
415 auto dof_tdir_arr = dofs_mfarrvec[lev][tdir]->const_array(mfi);
416#elif defined(WARPX_DIM_3D)
417 int tdir1 = (dir + 1) % 3;
418 int tdir2 = (dir + 2) % 3;
419 GpuArray<Array4<const int>, AMREX_SPACEDIM>
420 const dof_arrays {{ AMREX_D_DECL( dofs_mfarrvec[lev][dir]->const_array(mfi),
421 dofs_mfarrvec[lev][tdir1]->const_array(mfi),
422 dofs_mfarrvec[lev][tdir2]->const_array(mfi) ) }};
423#elif !defined(WARPX_DIM_1D_Z)
424 amrex::ignore_unused(dxi);
425 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble encountered unknown dimensionality");
426#endif
427
428 // Add the Maxwell's piece
429 if (thetaDt > 0.0) {
430 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
431 {
432 int ridx_l = dof_arr(i,j,k,0);
433 if (ridx_l < 0) { return; }
434
435 int icol = 0;
436 int ridx_g = dof_arr(i,j,k,1);
437
438 r_indices_g_ptr[ridx_l] = ridx_g;
439
440 {
441 int cidx_g_lhs = dof_arr(i,j,k,1);
442 amrex::Real val = 1.0;
443 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_lhs, val,
444 &c_indices_g_ptr[ridx_l*nnz_max],
445 &a_ij_ptr[ridx_l*nnz_max],
446 nnz_max, icol );
447 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
448 }
449
450#if defined(WARPX_DIM_1D_Z)
451 if (dir != 2) {
452 {
453 int cidx_g_rhs = dof_arr(i,j,k,1);
454 amrex::Real val = 2.0_rt*alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0);
455 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
456 &c_indices_g_ptr[ridx_l*nnz_max],
457 &a_ij_ptr[ridx_l*nnz_max],
458 nnz_max, icol );
459 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
460 }
461 {
462 int cidx_g_rhs = dof_arr(i-1,j,k,1);
463 amrex::Real val = -alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,1);
464 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
465 &c_indices_g_ptr[ridx_l*nnz_max],
466 &a_ij_ptr[ridx_l*nnz_max],
467 nnz_max, icol );
468 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
469 }
470 {
471 int cidx_g_rhs = dof_arr(i+1,j,k,1);
472 amrex::Real val = -alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,1);
473 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
474 &c_indices_g_ptr[ridx_l*nnz_max],
475 &a_ij_ptr[ridx_l*nnz_max],
476 nnz_max, icol );
477 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
478 }
479 }
480#elif defined(WARPX_DIM_XZ)
481 {
482 int cidx_g_rhs = dof_arr(i,j,k,1);
483 amrex::Real val = 0.0;
484 if (dir == 0) {
485 val += 2.0_rt*alpha * dxi[1]*dxi[1] * BC_mask_Edir_arr(i,j,k,0);
486 } else if (dir == 2) {
487 val += 2.0_rt*alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0);
488 } else if (dir == 1) {
489 val += 2.0_rt*alpha * ( dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0)
490 + dxi[1]*dxi[1] * BC_mask_Edir_arr(i,j,k,2) );
491 }
492 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
493 &c_indices_g_ptr[ridx_l*nnz_max],
494 &a_ij_ptr[ridx_l*nnz_max],
495 nnz_max, icol );
496 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
497 }
498 if ((dir == 0) || (dir == 2)) {
499 {
500 int cidx_g_rhs = (dir == 0 ? dof_arr(i,j-1,k,1) : dof_arr(i-1,j,k,1));
501 amrex::Real val = -alpha * dxi[dir==0?1:0] * dxi[dir==0?1:0] * BC_mask_Edir_arr(i,j,k,1);
502 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
503 &c_indices_g_ptr[ridx_l*nnz_max],
504 &a_ij_ptr[ridx_l*nnz_max],
505 nnz_max, icol );
506 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
507 }
508 {
509 int cidx_g_rhs = (dir == 0 ? dof_arr(i,j+1,k,1) : dof_arr(i+1,j,k,1));
510 amrex::Real val = -alpha * dxi[dir==0?1:0] * dxi[dir==0?1:0] * BC_mask_Edir_arr(i,j,k,1);
511 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
512 &c_indices_g_ptr[ridx_l*nnz_max],
513 &a_ij_ptr[ridx_l*nnz_max],
514 nnz_max, icol );
515 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
516 }
517 // The following four blocks are for
518 // dir = 0: d/dx(dEz/dz) at Ex(i,j) with Ex centered in x and nodal in z
519 // dir = 2: d/dz(dEx/dx) at Ez(i,j) with Ez centered in z and nodal in x
520 {
521 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i,j-1,k,1) : dof_tdir_arr(i-1,j,k,1));
522 amrex::Real val = alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
523 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
524 &c_indices_g_ptr[ridx_l*nnz_max],
525 &a_ij_ptr[ridx_l*nnz_max],
526 nnz_max, icol );
527 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
528 }
529 {
530 int cidx_g_rhs = dof_tdir_arr(i,j,k,1);
531 amrex::Real val = -alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
532 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
533 &c_indices_g_ptr[ridx_l*nnz_max],
534 &a_ij_ptr[ridx_l*nnz_max],
535 nnz_max, icol );
536 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
537 }
538 {
539 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j-1,k,1) : dof_tdir_arr(i-1,j+1,k,1));
540 amrex::Real val = -alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
541 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
542 &c_indices_g_ptr[ridx_l*nnz_max],
543 &a_ij_ptr[ridx_l*nnz_max],
544 nnz_max, icol );
545 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
546 }
547 {
548 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j,k,1) : dof_tdir_arr(i,j+1,k,1));
549 amrex::Real val = alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
550 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
551 &c_indices_g_ptr[ridx_l*nnz_max],
552 &a_ij_ptr[ridx_l*nnz_max],
553 nnz_max, icol );
554 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
555 }
556 } else if (dir==1) {
557 for (int jdir = 0; jdir <= 2; jdir+=2) {
558 {
559 int cidx_g_rhs = (jdir == 0 ? dof_arr(i-1,j,k,1) : dof_arr(i,j-1,k,1));
560 amrex::Real val = -alpha * dxi[jdir==0?0:1] * dxi[jdir==0?0:1] * BC_mask_Edir_arr(i,j,k,jdir+1);
561 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
562 &c_indices_g_ptr[ridx_l*nnz_max],
563 &a_ij_ptr[ridx_l*nnz_max],
564 nnz_max, icol );
565 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
566 }
567 {
568 int cidx_g_rhs = (jdir == 0 ? dof_arr(i+1,j,k,1) : dof_arr(i,j+1,k,1));
569 amrex::Real val = -alpha * dxi[jdir==0?0:1] * dxi[jdir==0?0:1] * BC_mask_Edir_arr(i,j,k,jdir+1);
570 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
571 &c_indices_g_ptr[ridx_l*nnz_max],
572 &a_ij_ptr[ridx_l*nnz_max],
573 nnz_max, icol );
574 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
575 }
576 }
577 }
578#elif defined(WARPX_DIM_3D)
579 amrex::IntVect dvec(AMREX_D_DECL(dir,tdir1,tdir2));
580 amrex::IntVect ic(AMREX_D_DECL(i,j,k));
581 {
582 int cidx_g_rhs = dof_arrays[0](ic,1);
583 amrex::Real val = 2.0_rt * alpha * ( dxi[dvec[1]]*dxi[dvec[1]] * BC_mask_Edir_arr(i,j,k,0)
584 + dxi[dvec[2]]*dxi[dvec[2]] * BC_mask_Edir_arr(i,j,k,3) );
585 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
586 &c_indices_g_ptr[ridx_l*nnz_max],
587 &a_ij_ptr[ridx_l*nnz_max],
588 nnz_max, icol );
589 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
590 }
591 for (int ctr = -1; ctr <= 1; ctr += 2) {
592 for (int tdir = 1; tdir <= 2; tdir++) {
593 auto iv = ic; iv[dvec[tdir]] += ctr;
594 int cidx_g_rhs = dof_arrays[0](iv,1);
595 const int comp_shift = (dvec[tdir] == tdir1) ? 0:3;
596 amrex::Real val = -alpha * dxi[dvec[tdir]]*dxi[dvec[tdir]] * BC_mask_Edir_arr(i,j,k,comp_shift+1);
597 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
598 &c_indices_g_ptr[ridx_l*nnz_max],
599 &a_ij_ptr[ridx_l*nnz_max],
600 nnz_max, icol );
601 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
602 }
603 }
604 for (int ctr_dir = -1; ctr_dir <= 0; ctr_dir++) {
605 for (int ctr_tdir = -1; ctr_tdir <= 0; ctr_tdir++) {
606 for (int tdir = 1; tdir <= 2; tdir++) {
607 auto iv = ic; iv[dvec[0]] += (ctr_dir+1); iv[dvec[tdir]] += ctr_tdir;
608 auto sign = std::copysign(1,ctr_dir) * std::copysign(1,ctr_tdir);
609 int cidx_g_rhs = dof_arrays[tdir](iv,1);
610 const int comp_shift = (dvec[tdir] == tdir1) ? 0:3;
611 amrex::Real val = amrex::Real(sign) * alpha * dxi[dvec[0]]*dxi[dvec[tdir]] * BC_mask_Edir_arr(i,j,k,comp_shift+2);
612 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
613 &c_indices_g_ptr[ridx_l*nnz_max],
614 &a_ij_ptr[ridx_l*nnz_max],
615 nnz_max, icol );
616 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
617 }
618 }
619 }
620#endif
621 num_nz_ptr[ridx_l] = icol;
622 });
623 }
624
625 // Add the mass matrix piece
626 // See Figure B.9 of JCP 491, 112383 (2023) for an illustrative diagram
627 // of the mass matrices (https://doi.org/10.1016/j.jcp.2023.112383).
628 //
629 // The coupling of Jx(i,j,k) to Ex(i+i0,j+j0,k+k0), where i0 ranges from
630 // -MM_width[0] to +MM_width[0], j0 ranges from -MM_width[1] to +MM_width[1], and k0 ranges
631 // from -MM_width[2] to +MM_width[2], is stored as components of m_bcoefs[dir=0] (mass matrices).
632 // Similarly for Jy/Ey (m_bcoefs[dir=1]) and Jz/Ez (m_bcoefs[dir=2]).
633 // The mapping to the components is given by:
634 // mm_comp = i0 + MM_width[0] + MM_comp[0]*(j0 + MM_width[1]) + (MM_comp[0] + MM_comp[1])*(k0 + MM_width[2])
635 if (m_include_mass_matrices) {
636
637 auto sigma_ii_arr = (*m_bcoefs)[lev][dir]->const_array(mfi);
638 amrex::GpuArray<int,3> MM_ncomp = {1,1,1};
639 amrex::GpuArray<int,3> MM_width = {0,0,0};
640 for (int space_dir=0; space_dir<AMREX_SPACEDIM; space_dir++) {
641 MM_ncomp[space_dir] = m_ops->GetMassMatricesPCnComp(dir,space_dir);
642 MM_width[space_dir] = (MM_ncomp[space_dir] - 1)/2;
643 }
644
645 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
646 {
647 int ridx_l = dof_arr(i,j,k,0);
648 if (ridx_l < 0) { return; }
649
650 const amrex::IntVect iv_base = IntVect(AMREX_D_DECL(i,j,k));
651 int icol = num_nz_ptr[ridx_l];
652
653 int mm_comp = 0;
654 for (int comp2 = 0; comp2 < MM_ncomp[2]; comp2++) {
655 [[maybe_unused]] const int kk0 = comp2 - MM_width[2];
656 for (int comp1 = 0; comp1 < MM_ncomp[1]; comp1++) {
657 [[maybe_unused]] const int jj0 = comp1 - MM_width[1];
658 for (int comp0 = 0; comp0 < MM_ncomp[0]; comp0++) {
659 const int ii0 = comp0 - MM_width[0]; // ..., -2, -1, 0, 1, 2, ...
660 const amrex::IntVect iv_shift = IntVect(AMREX_D_DECL(ii0, jj0, kk0));
661 if (full_bx.contains(iv_base + iv_shift)) {
662 int cidx_g_rhs = dof_arr(iv_base + iv_shift,1);
663 amrex::Real val = sigma_ii_arr(iv_base,mm_comp);
664 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
665 &c_indices_g_ptr[ridx_l*nnz_max],
666 &a_ij_ptr[ridx_l*nnz_max],
667 nnz_max, icol );
668 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
669 }
670 ++mm_comp;
671 }
672 }
673 }
674
675 num_nz_ptr[ridx_l] = icol;
676 });
677 }
678 Gpu::synchronize();
679 }
680 }
681
682 nnz_actual = std::max(nnz_actual, *(nnz_actual_d.copyToHost()));
683 }
684
685 amrex::ParallelDescriptor::ReduceIntMax(&nnz_actual, 1);
686 return (nnz_actual - nnz_max);
687}
688
689template <class T, class Ops>
690void MatrixPC<T,Ops>::Apply (T& a_x, const T& a_b)
691{
692 // Given a right-hand-side b, solve:
693 // A x = b
694 // where A is the linear operator, in this case, the curl-curl
695 // operator:
696 // A x = curl (alpha * curl (x) ) + beta * x
697
698 BL_PROFILE("MatrixPC::Apply()");
699 using namespace amrex;
700
701 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
702 IsDefined(),
703 "MatrixPC::Apply() called on undefined object" );
704 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
705 a_x.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
706 "MatrixPC::Apply() - a_x must be Efield_fp type");
707 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
708 a_b.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
709 "MatrixPC::Apply() - a_b must be Efield_fp type");
710
711 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Apply() - native matrix solvers not implemented. Use with external library, eg, PETSc.");
712
713}
714
715#endif
AMReX_Array.H
BL_PROFILE
#define BL_PROFILE(a)
AMREX_ALWAYS_ASSERT
#define AMREX_ALWAYS_ASSERT(EX)
AMREX_FORCE_INLINE
#define AMREX_FORCE_INLINE
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
AMREX_GPU_HOST_DEVICE
#define AMREX_GPU_HOST_DEVICE
pp
amrex::ParmParse pp
AMReX_MultiFab.H
AMReX_ParmParse.H
AMREX_D_DECL
#define AMREX_D_DECL(a, b, c)
AMReX_Vector.H
PhysicalSpecies::alpha
@ alpha
Definition SpeciesPhysicalProperties.H:18
TextMsg.H
WARPX_ABORT_WITH_MESSAGE
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition TextMsg.H:15
WARPX_ALWAYS_ASSERT_WITH_MESSAGE
#define WARPX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition TextMsg.H:13
MatrixPC::IsDefined
bool IsDefined() const override
Check if the nonlinear solver has been defined.
Definition MatrixPC.H:171
MatrixPC::Apply
void Apply(T &, const T &) override
Apply (solve) the preconditioner given a RHS.
Definition MatrixPC.H:690
MatrixPC::m_ops
Ops * m_ops
Definition MatrixPC.H:180
MatrixPC::m_verbose
bool m_verbose
Definition MatrixPC.H:178
MatrixPC::m_r_indices_g
amrex::Gpu::DeviceVector< int > m_r_indices_g
Definition MatrixPC.H:193
MatrixPC::MatrixPC
MatrixPC(const MatrixPC &)=delete
MatrixPC::m_c_indices_g
amrex::Gpu::DeviceVector< int > m_c_indices_g
Definition MatrixPC.H:195
MatrixPC::printParameters
void printParameters() const override
Print parameters.
Definition MatrixPC.H:212
MatrixPC::Update
void Update(const T &a_U) override
Update the preconditioner.
Definition MatrixPC.H:284
MatrixPC::Define
void Define(const T &, Ops *const) override
Define the preconditioner.
Definition MatrixPC.H:229
MatrixPC::m_num_nz
amrex::Gpu::DeviceVector< int > m_num_nz
Definition MatrixPC.H:194
MatrixPC::setName
void setName(const std::string &a_name) override
Set the name for screen output and parsing inputs.
Definition MatrixPC.H:173
MatrixPC::getPCMatrix
void getPCMatrix(amrex::Gpu::DeviceVector< int > &a_r_indices_g, amrex::Gpu::DeviceVector< int > &a_num_nz, amrex::Gpu::DeviceVector< int > &a_c_indices_g, amrex::Gpu::DeviceVector< RT > &a_a_ij, int &a_n, int &a_ncols_max) override
Get the sparse matrix form of the preconditioner.
Definition MatrixPC.H:128
MatrixPC::m_pc_diag_only
bool m_pc_diag_only
Definition MatrixPC.H:187
MatrixPC::m_bcoefs
const amrex::Vector< amrex::Array< amrex::MultiFab *, 3 > > * m_bcoefs
Definition MatrixPC.H:198
MatrixPC::m_ndofs_l
int m_ndofs_l
Definition MatrixPC.H:185
MatrixPC::m_num_realloc
int m_num_realloc
Definition MatrixPC.H:200
MatrixPC::MatrixPC
MatrixPC()=default
Default constructor.
MatrixPC::m_geom
amrex::Vector< amrex::Geometry > m_geom
Definition MatrixPC.H:183
MatrixPC::m_num_amr_levels
int m_num_amr_levels
Definition MatrixPC.H:182
MatrixPC::RT
typename T::value_type RT
Definition MatrixPC.H:80
MatrixPC::readParameters
void readParameters()
Read parameters.
Definition MatrixPC.H:221
MatrixPC::operator=
MatrixPC & operator=(const MatrixPC &)=delete
MatrixPC::m_ndofs_g
int m_ndofs_g
Definition MatrixPC.H:186
MatrixPC::m_include_mass_matrices
bool m_include_mass_matrices
Definition MatrixPC.H:189
MatrixPC::Assemble
int Assemble(const T &a_U)
Assemble the matrix.
Definition MatrixPC.H:317
MatrixPC::m_pc_mat_nnz
int m_pc_mat_nnz
Definition MatrixPC.H:188
MatrixPC::m_name
std::string m_name
Definition MatrixPC.H:191
MatrixPC::~MatrixPC
~MatrixPC() override=default
Default destructor.
MatrixPC::MatrixPC
MatrixPC(MatrixPC &&) noexcept=delete
MatrixPC::m_is_defined
bool m_is_defined
Definition MatrixPC.H:177
MatrixPC::m_a_ij
amrex::Gpu::DeviceVector< amrex::Real > m_a_ij
Definition MatrixPC.H:196
Preconditioner::m_dt
RT m_dt
Definition Preconditioner.H:117
Preconditioner::Preconditioner
Preconditioner()=default
Default constructor.
amrex::BoxND::contains
__host__ __device__ bool contains(const IntVectND< dim > &p) const noexcept
amrex::FabArray::const_array
Array4< typename FabArray< FAB >::value_type const > const_array(const MFIter &mfi) const noexcept
amrex::Gpu::Buffer::data
T const * data() const noexcept
amrex::Real
amrex_real Real
amrex::Gpu::DeviceVector
PODVector< T, ArenaAllocator< T > > DeviceVector
amrex::ParallelDescriptor::ReduceIntMax
void ReduceIntMax(int &)
MatrixPCUtils
Definition MatrixPC.H:24
MatrixPCUtils::insertOrAdd
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool insertOrAdd(const int a_cidx, const amrex::Real a_val, int *const a_cidxs, amrex::Real *const a_aij, const int a_nnz, int &a_ncol)
Definition MatrixPC.H:26
ablastr::constant::SI::c
constexpr auto c
vacuum speed of light [m/s]
Definition constant.H:153
ablastr::warn_manager::WMRecordWarning
void WMRecordWarning(const std::string &topic, const std::string &text, const WarnPriority &priority=WarnPriority::medium)
Helper function to abbreviate the call to WarnManager::GetInstance().RecordWarning (recording a warni...
Definition WarnManager.cpp:320
amrex::Gpu::Atomic::Max
__host__ __device__ AMREX_FORCE_INLINE T Max(T *const m, T const value) noexcept
amrex::Gpu::synchronize
void synchronize() noexcept
amrex::Gpu::copyAsync
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
amrex::Gpu::deviceToDevice
static constexpr DeviceToDevice deviceToDevice
amrex::Gpu::streamSynchronize
void streamSynchronize() noexcept
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
amrex::ParallelFor
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::Box
BoxND< 3 > Box
amrex::IntVect
IntVectND< 3 > IntVect
warpx::fields::FieldType::Efield_fp
@ Efield_fp
Definition Fields.H:96
Generated by doxygen 1.13.2