16
16
#include < AMReX_ParmParse.H>
17
17
18
18
#include < vector>
19
+ #include < istream>
20
+ #include < filesystem>
19
21
20
22
/* *
21
23
* \brief Newton method to solve nonlinear equation of form:
@@ -45,7 +47,8 @@ public:
45
47
void Solve ( Vec& a_U,
46
48
const Vec& a_b,
47
49
amrex::Real a_time,
48
- amrex::Real a_dt ) const override ;
50
+ amrex::Real a_dt,
51
+ int a_step) const override ;
49
52
50
53
void GetSolverParams ( amrex::Real& a_rtol,
51
54
amrex::Real& a_atol,
@@ -200,6 +203,34 @@ void NewtonSolver<Vec,Ops>::Define ( const Vec& a_U,
200
203
201
204
this ->m_is_defined = true ;
202
205
206
+ // Create diagnostic file and write header
207
+ if (!this ->m_diagnostic_file .empty () && amrex::ParallelDescriptor::IOProcessor ()) {
208
+
209
+ std::filesystem::path const diagnostic_path (this ->m_diagnostic_file );
210
+ std::filesystem::path const diagnostic_dir = diagnostic_path.parent_path ();
211
+ if (!diagnostic_dir.empty ()) {
212
+ std::filesystem::create_directories (diagnostic_dir);
213
+ }
214
+
215
+ std::ofstream diagnostic_file{this ->m_diagnostic_file , std::ofstream::out | std::ofstream::trunc };
216
+ int c = 0 ;
217
+ diagnostic_file << " #" ;
218
+ diagnostic_file << " [" << c++ << " ]step()" ;
219
+ diagnostic_file << " " ;
220
+ diagnostic_file << " [" << c++ << " ]time(s)" ;
221
+ diagnostic_file << " " ;
222
+ diagnostic_file << " [" << c++ << " ]iters" ;
223
+ diagnostic_file << " " ;
224
+ diagnostic_file << " [" << c++ << " ]norm_abs" ;
225
+ diagnostic_file << " " ;
226
+ diagnostic_file << " [" << c++ << " ]norm_rel" ;
227
+ diagnostic_file << " " ;
228
+ diagnostic_file << " [" << c++ << " ]gmres_iters" ;
229
+ diagnostic_file << " " ;
230
+ diagnostic_file << " [" << c++ << " ]gmres_last_res" ;
231
+ diagnostic_file << " \n " ;
232
+ diagnostic_file.close ();
233
+ }
203
234
}
204
235
205
236
template <class Vec , class Ops >
@@ -211,6 +242,7 @@ void NewtonSolver<Vec,Ops>::ParseParameters ()
211
242
pp_newton.query (" relative_tolerance" , m_rtol);
212
243
pp_newton.query (" max_iterations" , m_maxits);
213
244
pp_newton.query (" require_convergence" , m_require_convergence);
245
+ pp_newton.query (" diagnostic_file" , this ->m_diagnostic_file );
214
246
215
247
const amrex::ParmParse pp_gmres (" gmres" );
216
248
pp_gmres.query (" verbose_int" , m_gmres_verbose_int);
@@ -227,7 +259,8 @@ template <class Vec, class Ops>
227
259
void NewtonSolver<Vec,Ops>::Solve ( Vec& a_U,
228
260
const Vec& a_b,
229
261
amrex::Real a_time,
230
- amrex::Real a_dt ) const
262
+ amrex::Real a_dt,
263
+ int a_step) const
231
264
{
232
265
BL_PROFILE (" NewtonSolver::Solve()" );
233
266
WARPX_ALWAYS_ASSERT_WITH_MESSAGE (
@@ -248,6 +281,7 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec& a_U,
248
281
amrex::Real norm_rel = 0 .;
249
282
250
283
int iter;
284
+ int linear_solver_iters = 0 ;
251
285
for (iter = 0 ; iter < m_maxits;) {
252
286
253
287
// Compute residual: F(U) = U - b - R(U)
@@ -293,6 +327,7 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec& a_U,
293
327
// Solve linear system for Newton step [Jac]*dU = F
294
328
m_dU.zero ();
295
329
m_linear_solver->solve ( m_dU, m_F, m_gmres_rtol, m_gmres_atol );
330
+ linear_solver_iters += m_linear_solver->getNumIters ();
296
331
297
332
// Update solution
298
333
a_U -= m_dU;
@@ -321,6 +356,26 @@ void NewtonSolver<Vec,Ops>::Solve ( Vec& a_U,
321
356
}
322
357
}
323
358
359
+ if (!this ->m_diagnostic_file .empty () && amrex::ParallelDescriptor::IOProcessor ()) {
360
+ std::ofstream diagnostic_file{this ->m_diagnostic_file , std::ofstream::out | std::ofstream::app};
361
+ diagnostic_file << std::setprecision (14 );
362
+ diagnostic_file << a_step;
363
+ diagnostic_file << " " ;;
364
+ diagnostic_file << a_time;
365
+ diagnostic_file << " " ;;
366
+ diagnostic_file << iter;
367
+ diagnostic_file << " " ;;
368
+ diagnostic_file << norm_abs;
369
+ diagnostic_file << " " ;;
370
+ diagnostic_file << norm_rel;
371
+ diagnostic_file << " " ;;
372
+ diagnostic_file << linear_solver_iters;
373
+ diagnostic_file << " " ;;
374
+ diagnostic_file << m_linear_solver->getResidualNorm ();
375
+ diagnostic_file << " \n " ;
376
+ diagnostic_file.close ();
377
+ }
378
+
324
379
}
325
380
326
381
template <class Vec , class Ops >
0 commit comments