Skip to content

Commit 250479a

Browse files
eebassoeebasso
committed
Print variables in solve_bicgstab from development branch
1 parent a7afcba commit 250479a

File tree

1 file changed

+87
-0
lines changed

1 file changed

+87
-0
lines changed

‎Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,8 @@ public:
5252

5353
[[nodiscard]] int getNumIters () const noexcept { return iter; }
5454

55+
void PrintDots (const std::string& title, const char* names, const MF& x, const MF& y);
56+
5557
private:
5658

5759
MLLinOpT<MF>& Lp;
@@ -61,6 +63,7 @@ private:
6163
int verbose = 0;
6264
int maxiter = 100;
6365
IntVect nghost = IntVect(0);
66+
IntVect ngrow = IntVect(0);
6467
int iter = -1;
6568
};
6669

@@ -87,9 +90,26 @@ int
8790
MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
8891
{
8992
BL_PROFILE("MLCGSolver::bicgstab");
93+
std::string MLCGSolverTitle = "MLCGSolver_BiCGStab: ";
94+
std::string TitleHalf_Iter = "MLCGSolver_BiCGStab: Half Iter ";
95+
std::string TitleIteration = "MLCGSolver_BiCGStab: Iteration ";
96+
std::string TitlePost_Calc = "MLCGSolver_BiCGStab: Post_Calc ";
97+
98+
std::string Half_Iter = TitleHalf_Iter + std::to_string(0);
99+
std::string Iteration = TitleIteration + std::to_string(0);
100+
std::string Post_Calc = TitlePost_Calc + std::to_string(0);
101+
102+
amrex::Print() << "MLCGSolver_BiCGStab: nghost = " << nghost << "\n";
103+
amrex::Print() << "MLCGSolver_BiCGStab: ngrow = " << sol.nGrowVect() << "\n";
104+
105+
ngrow = sol.nGrowVect();
106+
verbose = 3;
90107

91108
const int ncomp = sol.nComp();
92109

110+
// auto& solArr = sol.arrays()[0];
111+
// solArr(3,0,0) = RT(0.0);
112+
93113
const BoxArray& ba = sol.boxArray();
94114
const DistributionMapping& dm = sol.DistributionMap();
95115
const auto& factory = sol.Factory();
@@ -98,6 +118,8 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
98118
MF sh(ba, dm, ncomp, sol.nGrowVect(), MFInfo(), factory);
99119
ph.setVal(RT(0.0));
100120
sh.setVal(RT(0.0));
121+
// PrintDots(Iteration,"Dots(ph,ph) initialized",ph,ph);
122+
// PrintDots(Iteration,"Dots(sh,sh) initialized",sh,sh);
101123

102124
MF sorig(ba, dm, ncomp, nghost, MFInfo(), factory);
103125
MF p (ba, dm, ncomp, nghost, MFInfo(), factory);
@@ -107,15 +129,23 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
107129
MF v (ba, dm, ncomp, nghost, MFInfo(), factory);
108130
MF t (ba, dm, ncomp, nghost, MFInfo(), factory);
109131

132+
PrintDots(Iteration,"Dots(sol,sol) before Lp.correctionResidual",sol,sol);
110133
Lp.correctionResidual(amrlev, mglev, r, sol, rhs, MLLinOpT<MF>::BCMode::Homogeneous);
134+
PrintDots(Iteration,"Dots(sol,sol) after Lp.correctionResidual",sol,sol);
111135

112136
// Then normalize
137+
PrintDots(Iteration,"Dots(r,r) before Lp.normalize",r,r);
113138
Lp.normalize(amrlev, mglev, r);
139+
PrintDots(Iteration,"Dots(r,r) after Lp.normalize",r,r);
114140

115141
sorig.LocalCopy(sol,0,0,ncomp,nghost);
142+
PrintDots(Iteration,"Dots(sorig,sorig)",sorig,sorig);
143+
116144
rh.LocalCopy (r ,0,0,ncomp,nghost);
145+
PrintDots(Iteration,"Dots(rh,rh)",rh,rh);
117146

118147
sol.setVal(RT(0.0));
148+
PrintDots(Iteration,"Dots(sol,sol) after sol.setVal(0)",sol,sol);
119149

120150
RT rnorm = norm_inf(r);
121151
const RT rnorm0 = rnorm;
@@ -141,26 +171,44 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
141171

142172
for (; iter <= maxiter; ++iter)
143173
{
174+
Half_Iter = TitleHalf_Iter + std::to_string(iter);
175+
Iteration = TitleIteration + std::to_string(iter);
176+
Post_Calc = TitlePost_Calc + std::to_string(iter-1);
177+
144178
const RT rho = dotxy(rh,r);
179+
PrintDots(Post_Calc,"Dots(rh,r)",rh,r);
180+
145181
if ( rho == 0 )
146182
{
147183
ret = 1; break;
148184
}
149185
if ( iter == 1 )
150186
{
151187
p.LocalCopy(r,0,0,ncomp,nghost);
188+
PrintDots(Post_Calc, "Dots(p,p)",p,p);
152189
}
153190
else
154191
{
155192
const RT beta = (rho/rho_1)*(alpha/omega);
193+
// amrex::Print() << Iteration << " beta: " << beta << '\n';
194+
PrintDots(Post_Calc, "Dots(p,p) before ",p,p);
156195
MF::Saxpy(p, -omega, v, 0, 0, ncomp, nghost); // p += -omega*v
196+
PrintDots(Post_Calc, "Dots(p,p) after Saxpy",p,p);
157197
MF::Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta*p
198+
PrintDots(Post_Calc, "Dots(p,p) after Xpay",p,p);
158199
}
159200
ph.LocalCopy(p,0,0,ncomp,nghost);
201+
PrintDots(Half_Iter,"Dots(ph,ph) before Lp.apply",ph,ph);
160202
Lp.apply(amrlev, mglev, v, ph, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
203+
PrintDots(Half_Iter,"Dots(ph,ph) after Lp.apply",ph,ph);
204+
205+
PrintDots(Half_Iter,"Dots(v,v) before Lp.normalize",v,v);
161206
Lp.normalize(amrlev, mglev, v);
207+
PrintDots(Half_Iter,"Dots(v,v) after Lp.normalize",v,v);
162208

163209
RT rhTv = dotxy(rh,v);
210+
PrintDots(Half_Iter,"Dots(rh,v)",rh,v);
211+
164212
if ( rhTv != RT(0.0) )
165213
{
166214
alpha = rho/rhTv;
@@ -169,8 +217,11 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
169217
{
170218
ret = 2; break;
171219
}
220+
// amrex::Print() << Half_Iter << " alpha: " << alpha << "\n";
172221
MF::Saxpy(sol, alpha, ph, 0, 0, ncomp, nghost); // sol += alpha * ph
173222
MF::LinComb(s, RT(1.0), r, 0, -alpha, v, 0, 0, ncomp, nghost); // s = r - alpha * v
223+
PrintDots(Half_Iter,"Dots(sol,sol)",sol,sol);
224+
PrintDots(Half_Iter,"Dots(s ,s )",s,s);
174225

175226
rnorm = norm_inf(s);
176227

@@ -185,8 +236,13 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
185236
if ( rnorm < eps_rel*rnorm0 || rnorm < eps_abs ) { break; }
186237

187238
sh.LocalCopy(s,0,0,ncomp,nghost);
239+
PrintDots(Iteration,"Dots(sh,sh) before Lp.apply",sh,sh);
188240
Lp.apply(amrlev, mglev, t, sh, MLLinOpT<MF>::BCMode::Homogeneous, MLLinOpT<MF>::StateMode::Correction);
241+
PrintDots(Iteration,"Dots(sh,sh) after Lp.apply",sh,sh);
242+
243+
PrintDots(Iteration,"Dots(t,t) before Lp.normalize",t,t);
189244
Lp.normalize(amrlev, mglev, t);
245+
PrintDots(Iteration,"Dots(t,t) after Lp.normalize",t,t);
190246
//
191247
// This is a little funky. I want to elide one of the reductions
192248
// in the following two dotxy()s. We do that by calculating the "local"
@@ -206,8 +262,12 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
206262
{
207263
ret = 3; break;
208264
}
265+
PrintDots(Iteration,"Dots(s,t)",t,s);
266+
209267
MF::Saxpy(sol, omega, sh, 0, 0, ncomp, nghost); // sol += omega * sh
210268
MF::LinComb(r, RT(1.0), s, 0, -omega, t, 0, 0, ncomp, nghost); // r = s - omega * t
269+
PrintDots(Iteration,"Dots(sol,sol)",sol,sol);
270+
PrintDots(Iteration,"Dots(r ,r )",r,r);
211271

212272
rnorm = norm_inf(r);
213273

@@ -244,14 +304,22 @@ MLCGSolverT<MF>::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs)
244304
ret = 8;
245305
}
246306

307+
std::string IterationFinal = "MLCGSolver_BiCGStab: Final: Iteration " + std::to_string(iter);
308+
247309
if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) )
248310
{
311+
// keep updated solution
312+
amrex::Print() << "MLCGSolver_BiCGstab: Keep new solution.\n";
249313
sol.LocalAdd(sorig, 0, 0, ncomp, nghost);
314+
PrintDots(IterationFinal,"Dots(sol,sol)",sol,sol);
250315
}
251316
else
252317
{
318+
// reset solution to original copy
319+
amrex::Print() << "MLCGSolver_BiCGstab: Reset solution.\n";
253320
sol.setVal(RT(0.0));
254321
sol.LocalAdd(sorig, 0, 0, ncomp, nghost);
322+
PrintDots(IterationFinal,"Dots(sol,sol)",sol,sol);
255323
}
256324

257325
return ret;
@@ -414,6 +482,25 @@ MLCGSolverT<MF>::norm_inf (const MF& res, bool local) -> RT
414482
return result;
415483
}
416484

485+
template <typename MF>
486+
void
487+
MLCGSolverT<MF>::PrintDots (const std::string& title, const char* name, const MF& x, const MF& y)
488+
{
489+
amrex::Print() << title << " " << name << " -> ";
490+
amrex::Print() << "dotxy: " ;
491+
amrex::Print() << dotxy(x,y);
492+
amrex::Print() << ", Dot(...,nghost): ";
493+
amrex::Print() << Dot(x,0,y,0,x.nComp(),nghost);
494+
amrex::Print() << ", Dot(...,ngrow): ";
495+
if (x.nGrowVect().allGE(ngrow) && y.nGrowVect().allGE(ngrow)) {
496+
amrex::Print() << Dot(x,0,y,0,x.nComp(),ngrow) << "\n";
497+
}
498+
else {
499+
amrex::Print() << "N/A" << "\n";
500+
}
501+
}
502+
503+
417504
using MLCGSolver = MLCGSolverT<MultiFab>;
418505

419506
}

0 commit comments

Comments
 (0)