44#include " DirectSystemSolve.hpp"
55#include < cmath>
66
7- namespace hpcReact
7+ namespace hpcReact
88{
99
1010namespace nonlinearSolvers
@@ -21,12 +21,13 @@ namespace internal
2121 * @return The Euclidean norm of the vector.
2222 */
2323template < int N >
24- HPCREACT_HOST_DEVICE
25- double norm ( double const (&r)[N])
24+ HPCREACT_HOST_DEVICE
25+ double norm ( double const (&r)[N] )
2626{
27- double sum = 0.0 ;
28- for ( int i = 0 ; i < N; ++i ) sum += r[i] * r[i];
29- return ::sqrt ( sum );
27+ double sum = 0.0 ;
28+ for ( int i = 0 ; i < N; ++i )
29+ sum += r[i] * r[i];
30+ return ::sqrt ( sum );
3031}
3132
3233/* *
@@ -37,10 +38,11 @@ double norm( double const (&r)[N])
3738 * @param dx The second vector, which will be added to the first vector.
3839 */
3940template < int N >
40- HPCREACT_HOST_DEVICE
41- void add (double (&x)[N], double const (&dx)[N])
41+ HPCREACT_HOST_DEVICE
42+ void add ( double (& x)[N], double const (&dx)[N] )
4243{
43- for (int i = 0 ; i < N; ++i) x[i] += dx[i];
44+ for ( int i = 0 ; i < N; ++i )
45+ x[i] += dx[i];
4446}
4547
4648/* *
@@ -51,16 +53,17 @@ void add(double (&x)[N], double const (&dx)[N])
5153 * @param value The scaling factor.
5254 */
5355template < int N >
54- HPCREACT_HOST_DEVICE
55- void scale ( double (&x)[N], double const value )
56+ HPCREACT_HOST_DEVICE
57+ void scale ( double (& x)[N], double const value )
5658{
57- for (int i = 0 ; i < N; ++i) x[i] *= value;
59+ for ( int i = 0 ; i < N; ++i )
60+ x[i] *= value;
5861}
5962
6063}
6164
6265namespace utils
63- {
66+ {
6467// LCOV_EXCL_START
6568
6669/* *
@@ -74,113 +77,115 @@ namespace utils
7477 */
7578template < int N >
7679HPCREACT_HOST_DEVICE
77- void print (double const (&J)[N][N], double const (&r)[N], double const (&dx)[N])
80+ void print ( double const (&J)[N][N], double const (&r)[N], double const (&dx)[N] )
7881{
79- printf (" =======================\n Jacobian matrix:\n =======================\n " );
80- printf (" RowID ColID Value\n " ); for ( int i = 0 ; i < N; ++i )
82+ printf ( " =======================\n Jacobian matrix:\n =======================\n " );
83+ printf ( " RowID ColID Value\n " ); for ( int i = 0 ; i < N; ++i )
84+ {
85+ for ( int j = 0 ; j < N; ++j )
8186 {
82- for ( int j = 0 ; j < N; ++j )
83- {
84- printf (" %10d%16d%27.16e\n " , i, j, J[i][j]);
85- }
87+ printf ( " %10d%16d%27.16e\n " , i, j, J[i][j] );
8688 }
89+ }
8790
88- printf (" \n =======================\n System right-hand side:\n =======================\n " );
89- printf (" RowID Value\n " );
91+ printf ( " \n =======================\n System right-hand side:\n =======================\n " );
92+ printf ( " RowID Value\n " );
9093
91- for ( int i = 0 ; i < N; ++i )
92- {
93- printf (" %10d%27.16e\n " , i, r[i]);
94- }
95-
96- printf (" \n =======================\n Delta update vector:\n =======================\n " );
97- printf (" RowID Value\n " );
98-
99- for ( int i = 0 ; i < N; ++i )
100- {
101- printf (" %10d%27.16e\n " , i, dx[i]);
102- }
94+ for ( int i = 0 ; i < N; ++i )
95+ {
96+ printf ( " %10d%27.16e\n " , i, r[i] );
97+ }
98+
99+ printf ( " \n =======================\n Delta update vector:\n =======================\n " );
100+ printf ( " RowID Value\n " );
101+
102+ for ( int i = 0 ; i < N; ++i )
103+ {
104+ printf ( " %10d%27.16e\n " , i, dx[i] );
105+ }
103106}
104107// LCOV_EXCL_STOP
105108
106109}
107110
108111template < int N,
109- typename REAL_TYPE,
110- typename ResidualFunc,
112+ typename REAL_TYPE,
113+ typename ResidualFunc,
111114 typename JacobianFunc >
112- HPCREACT_HOST_DEVICE
113- void newtonRaphson ( REAL_TYPE (&x)[N],
115+ HPCREACT_HOST_DEVICE
116+ void newtonRaphson ( REAL_TYPE (& x)[N],
114117 ResidualFunc computeResidual,
115118 JacobianFunc computeJacobian,
116119 int maxIters = 25,
117120 double tol = 1e-10 )
118121{
119- REAL_TYPE residual[N];
120- REAL_TYPE dx[N];
121- REAL_TYPE jacobian[N][N];
122+ REAL_TYPE residual[N];
123+ REAL_TYPE dx[N];
124+ REAL_TYPE jacobian[N][N];
122125
123- for ( int iter = 0 ; iter < maxIters; ++iter )
124- {
125- computeResidual ( x, residual );
126+ for ( int iter = 0 ; iter < maxIters; ++iter )
127+ {
128+ computeResidual ( x, residual );
126129
127- if ( internal::norm< N >( residual ) < tol ) return ;
130+ if ( internal::norm< N >( residual ) < tol )
131+ return ;
128132
129- computeJacobian ( x, jacobian );
133+ computeJacobian ( x, jacobian );
130134
131- solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
135+ solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
132136
133- internal::add< N >( x, dx );
134- }
137+ internal::add< N >( x, dx );
138+ }
135139}
136140
137141template < int N,
138- typename REAL_TYPE,
139- typename FUNCTION_TYPE >
140- HPCREACT_HOST_DEVICE
141- bool newtonRaphson ( REAL_TYPE (&x)[N],
142+ typename REAL_TYPE,
143+ typename FUNCTION_TYPE >
144+ HPCREACT_HOST_DEVICE
145+ bool newtonRaphson ( REAL_TYPE (& x)[N],
142146 FUNCTION_TYPE computeResidualAndJacobian,
143147 int maxIters = 12,
144148 double tol = 1e-10,
145149 bool const do_print = false )
146150{
147- REAL_TYPE residual[N]{};
148- REAL_TYPE dx[N]{};
149- REAL_TYPE jacobian[N][N]{};
150- bool isConverged = false ;
151+ REAL_TYPE residual[N]{};
152+ REAL_TYPE dx[N]{};
153+ REAL_TYPE jacobian[N][N]{};
154+ bool isConverged = false ;
155+
156+ for ( int iter = 0 ; iter < maxIters; ++iter )
157+ {
158+ computeResidualAndJacobian ( x, residual, jacobian );
159+
160+ double const norm = internal::norm< N >( residual );
161+
162+ printf ( " --Iter %d: Residual norm = %.12e\n " , iter, norm );
151163
152- for ( int iter = 0 ; iter < maxIters; ++iter )
164+ if ( norm < tol )
153165 {
154- computeResidualAndJacobian ( x, residual, jacobian );
155-
156- double const norm = internal::norm< N >( residual );
157-
158- printf ( " --Iter %d: Residual norm = %.12e\n " , iter, norm );
159-
160- if ( norm < tol ) {
161- printf ( " --Converged.\n " );
162- isConverged = true ;
163- break ;
164- }
165- internal::scale<N>( residual, -1.0 );
166-
167- if ( do_print )
168- {
169- utils::print ( jacobian, residual, dx ); // LCOV_EXCL_LINE
170- }
171-
172- solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
173- internal::add< N >( x, dx );
174-
166+ printf ( " --Converged.\n " );
167+ isConverged = true ;
168+ break ;
175169 }
176-
177- if ( !isConverged )
170+ internal::scale< N >( residual, -1.0 );
171+
172+ if ( do_print )
178173 {
179- printf ( " --Newton solver error: Max iterations reached without convergence. \n " );
174+ utils::print ( jacobian, residual, dx ); // LCOV_EXCL_LINE
180175 }
181176
182- return isConverged;
177+ solveNxN_pivoted< REAL_TYPE, N >( jacobian, residual, dx );
178+ internal::add< N >( x, dx );
179+
180+ }
181+
182+ if ( !isConverged )
183+ {
184+ printf ( " --Newton solver error: Max iterations reached without convergence.\n " );
185+ }
186+
187+ return isConverged;
183188}
184189
185190}
186- }
191+ }
0 commit comments