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 int m_pc_mass_matrix_width = 0;
190 bool m_pc_mass_matrix_include_ij = false;
191
192 std::string m_name = "noname";
193
194 amrex::Gpu::DeviceVector<int> m_r_indices_g;
195 amrex::Gpu::DeviceVector<int> m_num_nz;
196 amrex::Gpu::DeviceVector<int> m_c_indices_g;
197 amrex::Gpu::DeviceVector<amrex::Real> m_a_ij;
198
199 const amrex::Vector<amrex::Array<amrex::MultiFab*,3>>* m_bcoefs = nullptr;
200
201 int m_num_realloc = 0;
202
206 void readParameters();
207
208 private:
209
210};
211
212template <class T, class Ops>
213void MatrixPC<T,Ops>::printParameters() const
214{
215 using namespace amrex;
216 Print() << m_name << " verbose: " << (m_verbose?"true":"false") << "\n";
217 Print() << m_name << " pc_diagonal_only: " << (m_pc_diag_only?"true":"false") << "\n";
218 Print() << m_name << " pc_mass_matrix_width: " << m_pc_mass_matrix_width << "\n";
219 Print() << m_name << " pc_mass_matrix_include_ij: " << (m_pc_mass_matrix_include_ij?"true":"false") << "\n";
220}
221
222template <class T, class Ops>
223void MatrixPC<T,Ops>::readParameters()
224{
225 const amrex::ParmParse pp(m_name);
226 pp.query("verbose", m_verbose);
227 pp.query("pc_diagonal_only", m_pc_diag_only);
228 pp.query("pc_mass_matrix_width", m_pc_mass_matrix_width);
229 if(m_pc_mass_matrix_width>0) {
230 pp.query("pc_mass_matrix_include_ij", m_pc_mass_matrix_include_ij);
231 }
232}
233
234template <class T, class Ops>
235void MatrixPC<T,Ops>::Define ( const T& a_U,
236 Ops* const a_ops )
237{
238 BL_PROFILE("MatrixPC::Define()");
239 using namespace amrex;
240
241 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
242 !IsDefined(),
243 "MatrixPC::Define() called on defined object" );
244 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
245 (a_ops != nullptr),
246 "MatrixPC::Define(): a_ops is nullptr" );
247 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
248 a_U.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
249 "MatrixPC::Define() must be called with Efield_fp type");
250
251 m_ops = a_ops;
252 // read preconditioner parameters
253 readParameters();
254
255 // a_U is not needed
256 amrex::ignore_unused(a_U);
257
258 // Set number of AMR levels and create geometry, grids, and
259 // distribution mapping vectors.
260 m_num_amr_levels = m_ops->numAMRLevels();
261 if (m_num_amr_levels > 1) {
262 WARPX_ABORT_WITH_MESSAGE("MatrixPC::Define(): m_num_amr_levels > 1");
263 }
264 m_geom.resize(m_num_amr_levels);
265 for (int n = 0; n < m_num_amr_levels; n++) {
266 m_geom[n] = m_ops->GetGeometry(n);
267 }
268
269 m_ndofs_l = a_U.nDOF_local();
270 m_ndofs_g = a_U.nDOF_global();
271
272 auto n_rows = size_t(m_ndofs_l);
273 auto n_cols = size_t(m_pc_mat_nnz) * size_t(m_ndofs_l);
274
275 m_r_indices_g.resize(n_rows);
276 m_num_nz.resize(n_rows);
277 m_c_indices_g.resize(n_cols);
278 m_a_ij.resize(n_cols);
279
280
281 m_bcoefs = m_ops->GetMassMatricesCoeff();
282 if (m_bcoefs == nullptr) {
283 if (m_pc_mass_matrix_width >= 0) {
284 std::stringstream warning_message;
285 warning_message << "pc_mass_matrix_width is nonnegative but unable to get mass matrices from operator.\n";
286 warning_message << "setting m_pc_mass_matrix_width = -1 (no mass matrix in PC).\n";
287 ablastr::warn_manager::WMRecordWarning("MatrixPC", warning_message.str());
288 }
289 m_pc_mass_matrix_width = -1;
290 }
291
292 m_is_defined = true;
293}
294
295template <class T, class Ops>
296void MatrixPC<T,Ops>::Update (const T& a_U)
297{
298 BL_PROFILE("MatrixPC::Update()");
299 using namespace amrex;
300
301 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
302 IsDefined(),
303 "MatrixPC::Update() called on undefined object" );
304
305 while(1) {
306
307 auto nnz_diff = Assemble(a_U);
308 AMREX_ALWAYS_ASSERT(nnz_diff >= 0);
309 if (nnz_diff) {
310
311 m_pc_mat_nnz += nnz_diff;
312 m_num_realloc++;
313
314 } else {
315 break;
316 }
317 }
318
319 if (m_num_realloc > 1) {
320 std::stringstream warning_message;
321 warning_message << "Number of times arrays were reallocated due to new nonzero elements "
322 << "is greater than 1 (" << m_num_realloc <<"). This is unexpected.\n";
323 ablastr::warn_manager::WMRecordWarning("MatrixPC", warning_message.str());
324 }
325
326}
327
328template <class T, class Ops>
329int MatrixPC<T,Ops>::Assemble (const T& a_U)
330{
331 // Assemble the sparse matrix representation of the preconditioner
332 // A = curl (alpha * curl []) + M
333 // where M is the mass matrix. The following data is set in this function:
334 // - m_r_indices_g: integer array of size n with the global row indices
335 // - m_num_nz: integer array of size n with the number of non-zero elements
336 // in each row
337 // - m_c_indices: integer array of size n*ncmax with the global column indices
338 // of non-zero elements in each row (row-major)
339 // - m_a_ij: real-type array of size n*ncmax with the matrix element values
340 // (row-major format)
341 // where n is the local number of rows, and ncmax is the maximum number of non-zero
342 // elements per row.
343
344 BL_PROFILE("MatrixPC::Assemble()");
345 using namespace amrex;
346
347 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
348 IsDefined(),
349 "MatrixPC::Assemble() called on undefined object" );
350
351 // set the alpha coefficient for the curl-curl op
352 const RT thetaDt = m_ops->GetThetaForPC()*this->m_dt;
353 const RT alpha = (thetaDt*PhysConst::c) * (thetaDt*PhysConst::c);
354 if (m_verbose) {
355 amrex::Print() << "Updating " << m_name
356 << ": theta*dt = " << thetaDt << ", "
357 << " coefficients: "
358 << "alpha = " << alpha << "\n";
359 }
360
361 // Get DOF object from a_U
362 const auto& dofs_obj = a_U.getDOFsObject();
363 const auto& dofs_mfarrvec = dofs_obj->m_array;
364 AMREX_ALWAYS_ASSERT(m_ndofs_l == dofs_obj->m_nDoFs_l);
365 AMREX_ALWAYS_ASSERT(m_ndofs_g == dofs_obj->m_nDoFs_g);
366
367 m_r_indices_g.clear();
368 m_num_nz.clear();
369 m_c_indices_g.clear();
370 m_a_ij.clear();
371
372 auto n_rows = size_t(m_ndofs_l);
373 auto n_cols = size_t(m_pc_mat_nnz) * size_t(m_ndofs_l);
374
375 m_r_indices_g.resize(n_rows);
376 m_num_nz.resize(n_rows);
377 m_c_indices_g.resize(n_cols);
378 m_a_ij.resize(n_cols);
379
380 auto r_indices_g_ptr = m_r_indices_g.data();
381 auto num_nz_ptr = m_num_nz.data();
382 auto c_indices_g_ptr = m_c_indices_g.data();
383 auto a_ij_ptr = m_a_ij.data();
384
385 const auto nnz_max = m_pc_mat_nnz;
386 auto nnz_actual = nnz_max;
387
388 for (int lev = 0; lev < m_num_amr_levels; lev++) {
389
390 auto ncomp = dofs_mfarrvec[lev][0]->nComp();
391 AMREX_ALWAYS_ASSERT(ncomp == 2); // local, global
392
393 const auto& geom = m_geom[lev];
394 const auto dxi = geom.InvCellSizeArray();
395
396 Gpu::Buffer<int> nnz_actual_d({nnz_max});
397 auto* nnz_actual_ptr = nnz_actual_d.data();
398
399 for (int dir = 0; dir < 3; dir++) {
400
401 const amrex::MultiFab* BC_mask_Edir = m_ops->GetCurl2BCmask(lev,dir);
402
403 for (amrex::MFIter mfi(*dofs_mfarrvec[lev][dir]); mfi.isValid(); ++mfi) {
404
405 auto bx = mfi.tilebox();
406
407 // BC mask values for in-plane components of curl curl E operator for 2D
408 auto BC_mask_Edir_arr = BC_mask_Edir->const_array(mfi);
409
410 auto dof_arr = dofs_mfarrvec[lev][dir]->const_array(mfi);
411
412#if defined(WARPX_DIM_RCYLINDER)
413 amrex::ignore_unused(dxi, BC_mask_Edir_arr);
414 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RCYLINDER");
415#elif defined(WARPX_DIM_RSPHERE)
416 amrex::ignore_unused(dxi, BC_mask_Edir_arr);
417 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RSPHERE");
418#elif defined(WARPX_DIM_RZ)
419 amrex::ignore_unused(dxi, BC_mask_Edir_arr);
420 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RZ");
421#elif defined(WARPX_DIM_XZ)
422 int tdir = -1;
423 if (dir == 0) { tdir = 2; }
424 else if (dir == 2) { tdir = 0; }
425 else { tdir = 1; }
426 auto dof_tdir_arr = dofs_mfarrvec[lev][tdir]->const_array(mfi);
427#elif defined(WARPX_DIM_3D)
428 int tdir1 = (dir + 1) % 3;
429 int tdir2 = (dir + 2) % 3;
430 GpuArray<Array4<const int>, AMREX_SPACEDIM>
431 const dof_arrays {{ AMREX_D_DECL( dofs_mfarrvec[lev][dir]->const_array(mfi),
432 dofs_mfarrvec[lev][tdir1]->const_array(mfi),
433 dofs_mfarrvec[lev][tdir2]->const_array(mfi) ) }};
434#elif !defined(WARPX_DIM_1D_Z)
435 amrex::ignore_unused(dxi);
436 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble encountered unknown dimensionality");
437#endif
438
439 // Add the Maxwell's piece
440 if (thetaDt > 0.0) {
441 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
442 {
443 int ridx_l = dof_arr(i,j,k,0);
444 if (ridx_l < 0) { return; }
445
446 int icol = 0;
447 int ridx_g = dof_arr(i,j,k,1);
448
449 r_indices_g_ptr[ridx_l] = ridx_g;
450
451 {
452 int cidx_g_lhs = dof_arr(i,j,k,1);
453 amrex::Real val = 1.0;
454 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_lhs, val,
455 &c_indices_g_ptr[ridx_l*nnz_max],
456 &a_ij_ptr[ridx_l*nnz_max],
457 nnz_max, icol );
458 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
459 }
460
461#if defined(WARPX_DIM_1D_Z)
462 if (dir != 2) {
463 {
464 int cidx_g_rhs = dof_arr(i,j,k,1);
465 amrex::Real val = 2.0_rt*alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0);
466 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
467 &c_indices_g_ptr[ridx_l*nnz_max],
468 &a_ij_ptr[ridx_l*nnz_max],
469 nnz_max, icol );
470 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
471 }
472 {
473 int cidx_g_rhs = dof_arr(i-1,j,k,1);
474 amrex::Real val = -alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,1);
475 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
476 &c_indices_g_ptr[ridx_l*nnz_max],
477 &a_ij_ptr[ridx_l*nnz_max],
478 nnz_max, icol );
479 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
480 }
481 {
482 int cidx_g_rhs = dof_arr(i+1,j,k,1);
483 amrex::Real val = -alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,1);
484 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
485 &c_indices_g_ptr[ridx_l*nnz_max],
486 &a_ij_ptr[ridx_l*nnz_max],
487 nnz_max, icol );
488 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
489 }
490 }
491#elif defined(WARPX_DIM_XZ)
492 {
493 int cidx_g_rhs = dof_arr(i,j,k,1);
494 amrex::Real val = 0.0;
495 if (dir == 0) {
496 val += 2.0_rt*alpha * dxi[1]*dxi[1] * BC_mask_Edir_arr(i,j,k,0);
497 } else if (dir == 2) {
498 val += 2.0_rt*alpha * dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0);
499 } else {
500 val += 2.0_rt*alpha * ( dxi[0]*dxi[0] * BC_mask_Edir_arr(i,j,k,0)
501 + dxi[1]*dxi[1] * BC_mask_Edir_arr(i,j,k,2) );
502 }
503 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
504 &c_indices_g_ptr[ridx_l*nnz_max],
505 &a_ij_ptr[ridx_l*nnz_max],
506 nnz_max, icol );
507 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
508 }
509 if ((dir == 0) || (dir == 2)) {
510 {
511 int cidx_g_rhs = (dir == 0 ? dof_arr(i,j-1,k,1) : dof_arr(i-1,j,k,1));
512 amrex::Real val = -alpha * dxi[dir==0?1:0] * dxi[dir==0?1:0] * BC_mask_Edir_arr(i,j,k,1);
513 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
514 &c_indices_g_ptr[ridx_l*nnz_max],
515 &a_ij_ptr[ridx_l*nnz_max],
516 nnz_max, icol );
517 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
518 }
519 {
520 int cidx_g_rhs = (dir == 0 ? dof_arr(i,j+1,k,1) : dof_arr(i+1,j,k,1));
521 amrex::Real val = -alpha * dxi[dir==0?1:0] * dxi[dir==0?1:0] * BC_mask_Edir_arr(i,j,k,1);
522 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
523 &c_indices_g_ptr[ridx_l*nnz_max],
524 &a_ij_ptr[ridx_l*nnz_max],
525 nnz_max, icol );
526 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
527 }
528 // The following four blocks are for
529 // dir = 0: d/dx(dEz/dz) at Ex(i,j) with Ex centered in x and nodal in z
530 // dir = 2: d/dz(dEx/dx) at Ez(i,j) with Ez centered in z and nodal in x
531 {
532 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i,j-1,k,1) : dof_tdir_arr(i-1,j,k,1));
533 amrex::Real val = alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
534 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
535 &c_indices_g_ptr[ridx_l*nnz_max],
536 &a_ij_ptr[ridx_l*nnz_max],
537 nnz_max, icol );
538 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
539 }
540 {
541 int cidx_g_rhs = dof_tdir_arr(i,j,k,1);
542 amrex::Real val = -alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
543 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
544 &c_indices_g_ptr[ridx_l*nnz_max],
545 &a_ij_ptr[ridx_l*nnz_max],
546 nnz_max, icol );
547 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
548 }
549 {
550 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j-1,k,1) : dof_tdir_arr(i-1,j+1,k,1));
551 amrex::Real val = -alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
552 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
553 &c_indices_g_ptr[ridx_l*nnz_max],
554 &a_ij_ptr[ridx_l*nnz_max],
555 nnz_max, icol );
556 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
557 }
558 {
559 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j,k,1) : dof_tdir_arr(i,j+1,k,1));
560 amrex::Real val = alpha * dxi[0] * dxi[1] * BC_mask_Edir_arr(i,j,k,2);
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 } else {
568 for (int jdir = 0; jdir <= 2; jdir+=2) {
569 {
570 int cidx_g_rhs = (jdir == 0 ? dof_arr(i-1,j,k,1) : dof_arr(i,j-1,k,1));
571 amrex::Real val = -alpha * dxi[jdir==0?0:1] * dxi[jdir==0?0:1] * BC_mask_Edir_arr(i,j,k,jdir+1);
572 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
573 &c_indices_g_ptr[ridx_l*nnz_max],
574 &a_ij_ptr[ridx_l*nnz_max],
575 nnz_max, icol );
576 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
577 }
578 {
579 int cidx_g_rhs = (jdir == 0 ? dof_arr(i+1,j,k,1) : dof_arr(i,j+1,k,1));
580 amrex::Real val = -alpha * dxi[jdir==0?0:1] * dxi[jdir==0?0:1] * BC_mask_Edir_arr(i,j,k,jdir+1);
581 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
582 &c_indices_g_ptr[ridx_l*nnz_max],
583 &a_ij_ptr[ridx_l*nnz_max],
584 nnz_max, icol );
585 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
586 }
587 }
588 }
589#elif defined(WARPX_DIM_3D)
590 amrex::IntVect dvec(AMREX_D_DECL(dir,tdir1,tdir2));
591 amrex::IntVect ic(AMREX_D_DECL(i,j,k));
592 {
593 int cidx_g_rhs = dof_arrays[0](ic,1);
594 amrex::Real val = 2.0_rt * alpha * ( dxi[dvec[1]]*dxi[dvec[1]] * BC_mask_Edir_arr(i,j,k,0)
595 + dxi[dvec[2]]*dxi[dvec[2]] * BC_mask_Edir_arr(i,j,k,3) );
596 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
597 &c_indices_g_ptr[ridx_l*nnz_max],
598 &a_ij_ptr[ridx_l*nnz_max],
599 nnz_max, icol );
600 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
601 }
602 for (int ctr = -1; ctr <= 1; ctr += 2) {
603 for (int tdir = 1; tdir <= 2; tdir++) {
604 auto iv = ic; iv[dvec[tdir]] += ctr;
605 int cidx_g_rhs = dof_arrays[0](iv,1);
606 const int comp_shift = (dvec[tdir] == tdir1) ? 0:3;
607 amrex::Real val = -alpha * dxi[dvec[tdir]]*dxi[dvec[tdir]] * BC_mask_Edir_arr(i,j,k,comp_shift+1);
608 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
609 &c_indices_g_ptr[ridx_l*nnz_max],
610 &a_ij_ptr[ridx_l*nnz_max],
611 nnz_max, icol );
612 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
613 }
614 }
615 for (int ctr_dir = -1; ctr_dir <= 0; ctr_dir++) {
616 for (int ctr_tdir = -1; ctr_tdir <= 0; ctr_tdir++) {
617 for (int tdir = 1; tdir <= 2; tdir++) {
618 auto iv = ic; iv[dvec[0]] += (ctr_dir+1); iv[dvec[tdir]] += ctr_tdir;
619 auto sign = std::copysign(1,ctr_dir) * std::copysign(1,ctr_tdir);
620 int cidx_g_rhs = dof_arrays[tdir](iv,1);
621 const int comp_shift = (dvec[tdir] == tdir1) ? 0:3;
622 amrex::Real val = amrex::Real(sign) * alpha * dxi[dvec[0]]*dxi[dvec[tdir]] * BC_mask_Edir_arr(i,j,k,comp_shift+2);
623 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
624 &c_indices_g_ptr[ridx_l*nnz_max],
625 &a_ij_ptr[ridx_l*nnz_max],
626 nnz_max, icol );
627 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
628 }
629 }
630 }
631#endif
632 num_nz_ptr[ridx_l] = icol;
633 });
634 }
635
636 // Add the mass matrix piece
637 if (m_pc_mass_matrix_width >= 0) {
638
639 auto mm_width = m_pc_mass_matrix_width;
640 auto mm_incl_ij = m_pc_mass_matrix_include_ij;
641 auto sigma_ii_arr = (*m_bcoefs)[lev][dir]->const_array(mfi);
642
643 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
644 {
645 int ridx_l = dof_arr(i,j,k,0);
646 if (ridx_l < 0) { return; }
647
648 int icol = num_nz_ptr[ridx_l];
649
650 {
651 int cidx_g_lhs = dof_arr(i,j,k,1);
652 amrex::Real val = sigma_ii_arr(i,j,k,0);
653 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_lhs, val,
654 &c_indices_g_ptr[ridx_l*nnz_max],
655 &a_ij_ptr[ridx_l*nnz_max],
656 nnz_max, icol );
657 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
658 }
659
660 amrex::ignore_unused(mm_width); // will be used when adding off-diagonal elements
661 amrex::ignore_unused(mm_incl_ij); // will be used when adding off-diagonal elements
662
663 num_nz_ptr[ridx_l] = icol;
664 });
665 }
666 Gpu::synchronize();
667 }
668 }
669
670 nnz_actual = std::max(nnz_actual, *(nnz_actual_d.copyToHost()));
671 }
672
673 amrex::ParallelDescriptor::ReduceIntMax(&nnz_actual, 1);
674 return (nnz_actual - nnz_max);
675}
676
677template <class T, class Ops>
678void MatrixPC<T,Ops>::Apply (T& a_x, const T& a_b)
679{
680 // Given a right-hand-side b, solve:
681 // A x = b
682 // where A is the linear operator, in this case, the curl-curl
683 // operator:
684 // A x = curl (alpha * curl (x) ) + beta * x
685
686 BL_PROFILE("MatrixPC::Apply()");
687 using namespace amrex;
688
689 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
690 IsDefined(),
691 "MatrixPC::Apply() called on undefined object" );
692 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
693 a_x.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
694 "MatrixPC::Apply() - a_x must be Efield_fp type");
695 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
696 a_b.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
697 "MatrixPC::Apply() - a_b must be Efield_fp type");
698
699 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Apply() - native matrix solvers not implemented. Use with external library, eg, PETSc.");
700
701}
702
703#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:678
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:194
MatrixPC::MatrixPC
MatrixPC(const MatrixPC &)=delete
MatrixPC::m_c_indices_g
amrex::Gpu::DeviceVector< int > m_c_indices_g
Definition MatrixPC.H:196
MatrixPC::printParameters
void printParameters() const override
Print parameters.
Definition MatrixPC.H:213
MatrixPC::Update
void Update(const T &a_U) override
Update the preconditioner.
Definition MatrixPC.H:296
MatrixPC::Define
void Define(const T &, Ops *const) override
Define the preconditioner.
Definition MatrixPC.H:235
MatrixPC::m_num_nz
amrex::Gpu::DeviceVector< int > m_num_nz
Definition MatrixPC.H:195
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:199
MatrixPC::m_ndofs_l
int m_ndofs_l
Definition MatrixPC.H:185
MatrixPC::m_num_realloc
int m_num_realloc
Definition MatrixPC.H:201
MatrixPC::MatrixPC
MatrixPC()=default
Default constructor.
MatrixPC::m_geom
amrex::Vector< amrex::Geometry > m_geom
Definition MatrixPC.H:183
MatrixPC::m_pc_mass_matrix_include_ij
bool m_pc_mass_matrix_include_ij
Definition MatrixPC.H:190
MatrixPC::m_num_amr_levels
int m_num_amr_levels
Definition MatrixPC.H:182
MatrixPC::m_pc_mass_matrix_width
int m_pc_mass_matrix_width
Definition MatrixPC.H:189
MatrixPC::RT
typename T::value_type RT
Definition MatrixPC.H:80
MatrixPC::readParameters
void readParameters()
Read parameters.
Definition MatrixPC.H:223
MatrixPC::operator=
MatrixPC & operator=(const MatrixPC &)=delete
MatrixPC::m_ndofs_g
int m_ndofs_g
Definition MatrixPC.H:186
MatrixPC::Assemble
int Assemble(const T &a_U)
Assemble the matrix.
Definition MatrixPC.H:329
MatrixPC::m_pc_mat_nnz
int m_pc_mat_nnz
Definition MatrixPC.H:188
MatrixPC::m_name
std::string m_name
Definition MatrixPC.H:192
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:197
Preconditioner::m_dt
RT m_dt
Definition Preconditioner.H:117
Preconditioner::Preconditioner
Preconditioner()=default
Default constructor.
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::IntVect
IntVectND< 3 > IntVect
warpx::fields::FieldType::Efield_fp
@ Efield_fp
Definition Fields.H:96
Generated by doxygen 1.13.2