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 const auto& dofs_lhs_mfarrvec = dofs_obj->m_array_lhs;
365 AMREX_ALWAYS_ASSERT(m_ndofs_l == dofs_obj->m_nDoFs_l);
366 AMREX_ALWAYS_ASSERT(m_ndofs_g == dofs_obj->m_nDoFs_g);
367
368 m_r_indices_g.clear();
369 m_num_nz.clear();
370 m_c_indices_g.clear();
371 m_a_ij.clear();
372
373 auto n_rows = size_t(m_ndofs_l);
374 auto n_cols = size_t(m_pc_mat_nnz) * size_t(m_ndofs_l);
375
376 m_r_indices_g.resize(n_rows);
377 m_num_nz.resize(n_rows);
378 m_c_indices_g.resize(n_cols);
379 m_a_ij.resize(n_cols);
380
381 auto r_indices_g_ptr = m_r_indices_g.data();
382 auto num_nz_ptr = m_num_nz.data();
383 auto c_indices_g_ptr = m_c_indices_g.data();
384 auto a_ij_ptr = m_a_ij.data();
385
386 const auto ndofs_l = m_ndofs_l;
387 const auto ndofs_g = m_ndofs_g;
388 const auto nnz_max = m_pc_mat_nnz;
389 auto nnz_actual = nnz_max;
390
391 for (int lev = 0; lev < m_num_amr_levels; lev++) {
392
393 auto ncomp = dofs_mfarrvec[lev][0]->nComp();
394 AMREX_ALWAYS_ASSERT(ncomp == 2); // local, global
395
396 const auto& geom = m_geom[lev];
397 const auto dxi = geom.InvCellSizeArray();
398
399 Gpu::Buffer<int> nnz_actual_d({nnz_max});
400 auto* nnz_actual_ptr = nnz_actual_d.data();
401
402 for (int dir = 0; dir < 3; dir++) {
403
404 for (amrex::MFIter mfi(*dofs_mfarrvec[lev][dir]); mfi.isValid(); ++mfi) {
405
406 auto bx = mfi.tilebox();
407 auto dof_lhs_arr = dofs_lhs_mfarrvec[lev][dir]->const_array(mfi);
408
409#if defined(WARPX_DIM_1D_Z)
410 auto dof_arr = dofs_mfarrvec[lev][dir]->const_array(mfi);
411#elif defined(WARPX_DIM_RCYLINDER)
412 amrex::ignore_unused(dxi);
413 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RCYLINDER");
414#elif defined(WARPX_DIM_RSPHERE)
415 amrex::ignore_unused(dxi);
416 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RSPHERE");
417#elif defined(WARPX_DIM_RZ)
418 amrex::ignore_unused(dxi);
419 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Assemble() not yet implemented for WARPX_DIM_RZ");
420#elif defined(WARPX_DIM_XZ)
421 auto dof_arr = dofs_mfarrvec[lev][dir]->const_array(mfi);
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#else
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 icol = 0;
444 int ridx_l = dof_lhs_arr(i,j,k,0);
445 int ridx_g = dof_lhs_arr(i,j,k,1);
446 AMREX_ALWAYS_ASSERT((ridx_l >= 0) && (ridx_l < ndofs_l));
447 AMREX_ALWAYS_ASSERT((ridx_g >= 0) && (ridx_g < ndofs_g));
448
449 r_indices_g_ptr[ridx_l] = ridx_g;
450
451 {
452 int cidx_g_lhs = dof_lhs_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];
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];
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];
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];
497 } else if (dir == 2) {
498 val += 2.0_rt*alpha * dxi[0]*dxi[0];
499 } else {
500 val += 2.0_rt*alpha *(dxi[0]*dxi[0] + dxi[1]*dxi[1]);
501 }
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 if ((dir == 0) || (dir == 2)) {
509 {
510 int cidx_g_rhs = (dir == 0 ? dof_arr(i,j-1,k,1) : dof_arr(i-1,j,k,1));
511 amrex::Real val = -alpha * dxi[dir==0?1:0] * dxi[dir==0?1:0];
512 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
513 &c_indices_g_ptr[ridx_l*nnz_max],
514 &a_ij_ptr[ridx_l*nnz_max],
515 nnz_max, icol );
516 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
517 }
518 {
519 int cidx_g_rhs = (dir == 0 ? dof_arr(i,j+1,k,1) : dof_arr(i+1,j,k,1));
520 amrex::Real val = -alpha * dxi[dir==0?1:0] * dxi[dir==0?1:0];
521 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
522 &c_indices_g_ptr[ridx_l*nnz_max],
523 &a_ij_ptr[ridx_l*nnz_max],
524 nnz_max, icol );
525 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
526 }
527 // The following four blocks are for
528 // dir = 0: d/dx(dEz/dz) at Ex(i,j) with Ex centered in x and nodal in z
529 // dir = 2: d/dz(dEx/dx) at Ez(i,j) with Ez centered in z and nodal in x
530 {
531 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i,j-1,k,1) : dof_tdir_arr(i-1,j,k,1));
532 amrex::Real val = alpha * dxi[0] * dxi[1];
533 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
534 &c_indices_g_ptr[ridx_l*nnz_max],
535 &a_ij_ptr[ridx_l*nnz_max],
536 nnz_max, icol );
537 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
538 }
539 {
540 int cidx_g_rhs = dof_tdir_arr(i,j,k,1);
541 amrex::Real val = -alpha * dxi[0] * dxi[1];
542 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
543 &c_indices_g_ptr[ridx_l*nnz_max],
544 &a_ij_ptr[ridx_l*nnz_max],
545 nnz_max, icol );
546 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
547 }
548 {
549 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j-1,k,1) : dof_tdir_arr(i-1,j+1,k,1));
550 amrex::Real val = -alpha * dxi[0] * dxi[1];
551 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
552 &c_indices_g_ptr[ridx_l*nnz_max],
553 &a_ij_ptr[ridx_l*nnz_max],
554 nnz_max, icol );
555 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
556 }
557 {
558 int cidx_g_rhs = (dir == 0 ? dof_tdir_arr(i+1,j,k,1) : dof_tdir_arr(i,j+1,k,1));
559 amrex::Real val = alpha * dxi[0] * dxi[1];
560 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
561 &c_indices_g_ptr[ridx_l*nnz_max],
562 &a_ij_ptr[ridx_l*nnz_max],
563 nnz_max, icol );
564 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
565 }
566 } else {
567 for (int jdir = 0; jdir <= 2; jdir+=2) {
568 {
569 int cidx_g_rhs = (jdir == 0 ? dof_arr(i-1,j,k,1) : dof_arr(i,j-1,k,1));
570 amrex::Real val = -alpha * dxi[jdir==0?0:1] * dxi[jdir==0?0:1];
571 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
572 &c_indices_g_ptr[ridx_l*nnz_max],
573 &a_ij_ptr[ridx_l*nnz_max],
574 nnz_max, icol );
575 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
576 }
577 {
578 int cidx_g_rhs = (jdir == 0 ? dof_arr(i+1,j,k,1) : dof_arr(i,j+1,k,1));
579 amrex::Real val = -alpha * dxi[jdir==0?0:1] * dxi[jdir==0?0:1];
580 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
581 &c_indices_g_ptr[ridx_l*nnz_max],
582 &a_ij_ptr[ridx_l*nnz_max],
583 nnz_max, icol );
584 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
585 }
586 }
587 }
588#elif defined(WARPX_DIM_3D)
589 amrex::IntVect dvec(AMREX_D_DECL(dir,tdir1,tdir2));
590 amrex::IntVect ic(AMREX_D_DECL(i,j,k));
591 {
592 int cidx_g_rhs = dof_arrays[0](ic,1);
593 amrex::Real val = 2.0_rt * alpha * (dxi[dvec[1]]*dxi[dvec[1]] + dxi[dvec[2]]*dxi[dvec[2]]);
594 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
595 &c_indices_g_ptr[ridx_l*nnz_max],
596 &a_ij_ptr[ridx_l*nnz_max],
597 nnz_max, icol );
598 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
599 }
600 for (int ctr = -1; ctr <= 1; ctr += 2) {
601 for (int tdir = 1; tdir <= 2; tdir++) {
602 auto iv = ic; iv[dvec[tdir]] += ctr;
603 int cidx_g_rhs = dof_arrays[0](iv,1);
604 amrex::Real val = -alpha * dxi[dvec[tdir]]*dxi[dvec[tdir]];
605 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
606 &c_indices_g_ptr[ridx_l*nnz_max],
607 &a_ij_ptr[ridx_l*nnz_max],
608 nnz_max, icol );
609 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
610 }
611 }
612 for (int ctr_dir = -1; ctr_dir <= 0; ctr_dir++) {
613 for (int ctr_tdir = -1; ctr_tdir <= 0; ctr_tdir++) {
614 for (int tdir = 1; tdir <= 2; tdir++) {
615 auto iv = ic; iv[dvec[0]] += (ctr_dir+1); iv[dvec[tdir]] += ctr_tdir;
616 auto sign = std::copysign(1,ctr_dir) * std::copysign(1,ctr_tdir);
617 int cidx_g_rhs = dof_arrays[tdir](iv,1);
618 amrex::Real val = amrex::Real(sign) * alpha * dxi[dvec[0]]*dxi[dvec[tdir]];
619 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_rhs, val,
620 &c_indices_g_ptr[ridx_l*nnz_max],
621 &a_ij_ptr[ridx_l*nnz_max],
622 nnz_max, icol );
623 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
624 }
625 }
626 }
627#endif
628 num_nz_ptr[ridx_l] = icol;
629 });
630 }
631
632 // Add the mass matrix piece
633 if (m_pc_mass_matrix_width >= 0) {
634
635 auto mm_width = m_pc_mass_matrix_width;
636 auto mm_incl_ij = m_pc_mass_matrix_include_ij;
637 auto sigma_ii_arr = (*m_bcoefs)[lev][dir]->const_array(mfi);
638
639 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
640 {
641 int ridx_l = dof_lhs_arr(i,j,k,0);
642 AMREX_ALWAYS_ASSERT((ridx_l >= 0) && (ridx_l < ndofs_l));
643 int icol = num_nz_ptr[ridx_l];
644
645 {
646 int cidx_g_lhs = dof_lhs_arr(i,j,k,1);
647 amrex::Real val = sigma_ii_arr(i,j,k,0) - 1.0;
648 auto flag = MatrixPCUtils::insertOrAdd( cidx_g_lhs, val,
649 &c_indices_g_ptr[ridx_l*nnz_max],
650 &a_ij_ptr[ridx_l*nnz_max],
651 nnz_max, icol );
652 if (!flag) { Gpu::Atomic::Max(nnz_actual_ptr, icol); }
653 }
654
655 amrex::ignore_unused(mm_width); // will be used when adding off-diagonal elements
656 amrex::ignore_unused(mm_incl_ij); // will be used when adding off-diagonal elements
657
658 num_nz_ptr[ridx_l] = icol;
659 });
660 }
661 Gpu::synchronize();
662 }
663 }
664
665 nnz_actual = std::max(nnz_actual, *(nnz_actual_d.copyToHost()));
666 }
667
668 amrex::ParallelDescriptor::ReduceIntMax(&nnz_actual, 1);
669 return (nnz_actual - nnz_max);
670}
671
672template <class T, class Ops>
673void MatrixPC<T,Ops>::Apply (T& a_x, const T& a_b)
674{
675 // Given a right-hand-side b, solve:
676 // A x = b
677 // where A is the linear operator, in this case, the curl-curl
678 // operator:
679 // A x = curl (alpha * curl (x) ) + beta * x
680
681 BL_PROFILE("MatrixPC::Apply()");
682 using namespace amrex;
683
684 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
685 IsDefined(),
686 "MatrixPC::Apply() called on undefined object" );
687 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
688 a_x.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
689 "MatrixPC::Apply() - a_x must be Efield_fp type");
690 WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
691 a_b.getArrayVecType()==warpx::fields::FieldType::Efield_fp,
692 "MatrixPC::Apply() - a_b must be Efield_fp type");
693
694 WARPX_ABORT_WITH_MESSAGE("MatrixPC<T,Ops>::Apply() - native matrix solvers not implemented. Use with external library, eg, PETSc.");
695
696}
697
698#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:673
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:125
Preconditioner::Preconditioner
Preconditioner()=default
Default constructor.
amrex::Gpu::Buffer::data
T const * data() const noexcept
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:149
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::Gpu::DeviceVector
PODVector< T, ArenaAllocator< T > > DeviceVector
amrex::ParallelDescriptor::ReduceIntMax
void ReduceIntMax(int &rvar)
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:94
Generated by doxygen 1.13.2