Peano
Loading...
Searching...
No Matches
Newton.cpp
Go to the documentation of this file.
1/* TwoPunctures: File "Newton.c"*/
2
3#include <stdio.h>
4#include <stdlib.h>
5#include <string.h>
6#include <math.h>
7#include <ctype.h>
8#include <gsl/gsl_vector.h>
9#include <gsl/gsl_linalg.h>
10#include "libtwopunctures/TwoPunctures.h"
11
12namespace TP {
13using namespace Utilities;
14
15const int StencilSize = 19;
16const int N_PlaneRelax = 1;
17const int NRELAX = 200;
18const int Step_Relax = 1;
19
20/* --------------------------------------------------------------------------*/
21double
22TwoPunctures::norm_inf (double const * TP_RESTRICT const F,
23 int const ntotal)
24{
25 double dmax = -1;
26 {
27 double dmax1 = -1;
28 for (int j = 0; j < ntotal; j++)
29 if (fabs (F[j]) > dmax1)
30 dmax1 = fabs (F[j]);
31 if (dmax1 > dmax)
32 dmax = dmax1;
33 }
34 return dmax;
35}
36/* --------------------------------------------------------------------------*/
37void
38TwoPunctures::resid (double * TP_RESTRICT const res,
39 int const ntotal,
40 double const * TP_RESTRICT const dv,
41 double const * TP_RESTRICT const rhs,
42 int const * TP_RESTRICT const ncols,
43 int const * TP_RESTRICT const * TP_RESTRICT const cols,
44 double const * TP_RESTRICT const * TP_RESTRICT const JFD)
45{
46 for (int i = 0; i < ntotal; i++)
47 {
48 double JFDdv_i = 0;
49 for (int m = 0; m < ncols[i]; m++)
50 JFDdv_i += JFD[i][m] * dv[cols[i][m]];
51 res[i] = rhs[i] - JFDdv_i;
52 }
53}
54
55/* -------------------------------------------------------------------------*/
56void
57TwoPunctures::LineRelax_al (double * TP_RESTRICT const dv,
58 int const j, int const k, int const nvar,
59 int const n1, int const n2, int const n3,
60 double const * TP_RESTRICT const rhs,
61 int const * TP_RESTRICT const ncols,
62 int const * TP_RESTRICT const * TP_RESTRICT const cols,
63 double const * TP_RESTRICT const * TP_RESTRICT const JFD)
64{
65 int i, m, Ic, Ip, Im, col, ivar;
66
67 gsl_vector *diag = gsl_vector_alloc(n1);
68 gsl_vector *e = gsl_vector_alloc(n1-1); /* above diagonal */
69 gsl_vector *f = gsl_vector_alloc(n1-1); /* below diagonal */
70 gsl_vector *b = gsl_vector_alloc(n1); /* rhs */
71 gsl_vector *x = gsl_vector_alloc(n1); /* solution vector */
72
73 for (ivar = 0; ivar < nvar; ivar++)
74 {
75 gsl_vector_set_zero(diag);
76 gsl_vector_set_zero(e);
77 gsl_vector_set_zero(f);
78 for (i = 0; i < n1; i++)
79 {
80 Ip = Index (ivar, i + 1, j, k, nvar, n1, n2, n3);
81 Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
82 Im = Index (ivar, i - 1, j, k, nvar, n1, n2, n3);
83 gsl_vector_set(b,i,rhs[Ic]);
84 for (m = 0; m < ncols[Ic]; m++)
85 {
86 col = cols[Ic][m];
87 if (col != Ip && col != Ic && col != Im)
88 *gsl_vector_ptr(b, i) -= JFD[Ic][m] * dv[col];
89 else
90 {
91 if (col == Im && i > 0)
92 gsl_vector_set(f,i-1,JFD[Ic][m]);
93 if (col == Ic)
94 gsl_vector_set(diag,i,JFD[Ic][m]);
95 if (col == Ip && i < n1-1)
96 gsl_vector_set(e,i,JFD[Ic][m]);
97 }
98 }
99 }
100 gsl_linalg_solve_tridiag(diag, e, f, b, x);
101 for (i = 0; i < n1; i++)
102 {
103 Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
104 dv[Ic] = gsl_vector_get(x, i);
105 }
106 }
107
108 gsl_vector_free(diag);
109 gsl_vector_free(e);
110 gsl_vector_free(f);
111 gsl_vector_free(b);
112 gsl_vector_free(x);
113}
114
115/* --------------------------------------------------------------------------*/
116void
117TwoPunctures::LineRelax_be (double * TP_RESTRICT const dv,
118 int const i, int const k, int const nvar,
119 int const n1, int const n2, int const n3,
120 double const * TP_RESTRICT const rhs,
121 int const * TP_RESTRICT const ncols,
122 int const * TP_RESTRICT const * TP_RESTRICT const cols,
123 double const * TP_RESTRICT const * TP_RESTRICT const JFD)
124{
125 int j, m, Ic, Ip, Im, col, ivar;
126
127 gsl_vector *diag = gsl_vector_alloc(n2);
128 gsl_vector *e = gsl_vector_alloc(n2-1); /* above diagonal */
129 gsl_vector *f = gsl_vector_alloc(n2-1); /* below diagonal */
130 gsl_vector *b = gsl_vector_alloc(n2); /* rhs */
131 gsl_vector *x = gsl_vector_alloc(n2); /* solution vector */
132
133 for (ivar = 0; ivar < nvar; ivar++)
134 {
135 gsl_vector_set_zero(diag);
136 gsl_vector_set_zero(e);
137 gsl_vector_set_zero(f);
138 for (j = 0; j < n2; j++)
139 {
140 Ip = Index (ivar, i, j + 1, k, nvar, n1, n2, n3);
141 Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
142 Im = Index (ivar, i, j - 1, k, nvar, n1, n2, n3);
143 gsl_vector_set(b,j,rhs[Ic]);
144 for (m = 0; m < ncols[Ic]; m++)
145 {
146 col = cols[Ic][m];
147 if (col != Ip && col != Ic && col != Im)
148 *gsl_vector_ptr(b, j) -= JFD[Ic][m] * dv[col];
149 else
150 {
151 if (col == Im && j > 0)
152 gsl_vector_set(f,j-1,JFD[Ic][m]);
153 if (col == Ic)
154 gsl_vector_set(diag,j,JFD[Ic][m]);
155 if (col == Ip && j < n2-1)
156 gsl_vector_set(e,j,JFD[Ic][m]);
157 }
158 }
159 }
160 gsl_linalg_solve_tridiag(diag, e, f, b, x);
161 for (j = 0; j < n2; j++)
162 {
163 Ic = Index (ivar, i, j, k, nvar, n1, n2, n3);
164 dv[Ic] = gsl_vector_get(x, j);
165 }
166 }
167 gsl_vector_free(diag);
168 gsl_vector_free(e);
169 gsl_vector_free(f);
170 gsl_vector_free(b);
171 gsl_vector_free(x);
172}
173
174/* --------------------------------------------------------------------------*/
175void
176TwoPunctures::relax (double * TP_RESTRICT const dv,
177 int const nvar, int const n1, int const n2, int const n3,
178 double const * TP_RESTRICT const rhs,
179 int const * TP_RESTRICT const ncols,
180 int const * TP_RESTRICT const * TP_RESTRICT const cols,
181 double const * TP_RESTRICT const * TP_RESTRICT const JFD)
182{
183 int i, j, k, n;
184
185 for (k = 0; k < n3; k = k + 2)
186 {
187 for (n = 0; n < N_PlaneRelax; n++)
188 {
189
190 for (i = 2; i < n1; i = i + 2)
191 LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
192
193 for (i = 1; i < n1; i = i + 2)
194 LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
195
196 for (j = 1; j < n2; j = j + 2)
197 LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
198
199 for (j = 0; j < n2; j = j + 2)
200 LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
201 }
202 }
203 for (k = 1; k < n3; k = k + 2)
204 {
205 for (n = 0; n < N_PlaneRelax; n++)
206 {
207
208 for (i = 0; i < n1; i = i + 2)
209 LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
210
211 for (i = 1; i < n1; i = i + 2)
212 LineRelax_be (dv, i, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
213
214 for (j = 1; j < n2; j = j + 2)
215 LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
216
217 for (j = 0; j < n2; j = j + 2)
218 LineRelax_al (dv, j, k, nvar, n1, n2, n3, rhs, ncols, cols, JFD);
219 }
220 }
221}
222
223/* --------------------------------------------------------------------------*/
224void
225TwoPunctures::TestRelax (int nvar, int n1, int n2, int n3, derivs v,
226 double *dv)
227{
228
229 int ntotal = n1 * n2 * n3 * nvar, **cols, *ncols, maxcol =
230 StencilSize * nvar, j;
231 double *F, *res, **JFD;
232 derivs u;
233
234 F = dvector (0, ntotal - 1);
235 res = dvector (0, ntotal - 1);
236 allocate_derivs (&u, ntotal);
237
238 JFD = dmatrix (0, ntotal - 1, 0, maxcol - 1);
239 cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
240 ncols = ivector (0, ntotal - 1);
241
242 F_of_v (nvar, n1, n2, n3, v, F, u);
243
244 SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
245
246 for (j = 0; j < ntotal; j++)
247 dv[j] = 0;
248 resid (res, ntotal, dv, F, ncols, cols, JFD);
249 printf ("Before: |F|=%20.15e\n", (double) norm1 (res, ntotal));
250 fflush(stdout);
251 for (j = 0; j < NRELAX; j++)
252 {
253 relax (dv, nvar, n1, n2, n3, F, ncols, cols, JFD); /* solves JFD*sh = s*/
254 if (j % Step_Relax == 0)
255 {
256 resid (res, ntotal, dv, F, ncols, cols, JFD);
257 printf ("j=%d\t |F|=%20.15e\n", j, (double) norm1 (res, ntotal));
258 fflush(stdout);
259 }
260 }
261
262 resid (res, ntotal, dv, F, ncols, cols, JFD);
263 printf ("After: |F|=%20.15e\n", (double) norm1 (res, ntotal));
264 fflush(stdout);
265
266 free_dvector (F, 0, ntotal - 1);
267 free_dvector (res, 0, ntotal - 1);
268 free_derivs (&u, ntotal);
269
270 free_dmatrix (JFD, 0, ntotal - 1, 0, maxcol - 1);
271 free_imatrix (cols, 0, ntotal - 1, 0, maxcol - 1);
272 free_ivector (ncols, 0, ntotal - 1);
273}
274
275/* --------------------------------------------------------------------------*/
276int
277TwoPunctures::bicgstab (int const nvar, int const n1, int const n2, int const n3,
278 derivs v, derivs dv,
279 int const output, int const itmax, double const tol,
280 double * TP_RESTRICT const normres)
281{
282
283 int ntotal = n1 * n2 * n3 * nvar, ii;
284 double alpha = 0, beta = 0;
285 double rho = 0, rho1 = 1, rhotol = 1e-50;
286 double omega = 0, omegatol = 1e-50;
287 double *p, *rt, *s, *t, *r, *vv;
288 double **JFD;
289 int **cols, *ncols, maxcol = StencilSize * nvar;
290 double *F;
291 derivs u, ph, sh;
292
293 F = dvector (0, ntotal - 1);
294 allocate_derivs (&u, ntotal);
295
296 JFD = dmatrix (0, ntotal - 1, 0, maxcol - 1);
297 cols = imatrix (0, ntotal - 1, 0, maxcol - 1);
298 ncols = ivector (0, ntotal - 1);
299
300 F_of_v (nvar, n1, n2, n3, v, F, u);
301 SetMatrix_JFD (nvar, n1, n2, n3, u, ncols, cols, JFD);
302
303 /* temporary storage */
304 r = dvector (0, ntotal - 1);
305 p = dvector (0, ntotal - 1);
306 allocate_derivs (&ph, ntotal);
307/* ph = dvector(0, ntotal-1);*/
308 rt = dvector (0, ntotal - 1);
309 s = dvector (0, ntotal - 1);
310 allocate_derivs (&sh, ntotal);
311/* sh = dvector(0, ntotal-1);*/
312 t = dvector (0, ntotal - 1);
313 vv = dvector (0, ntotal - 1);
314
315 /* check */
316 if (output == 1) {
317 printf ("bicgstab: itmax %d, tol %e\n", itmax, (double)tol);
318 fflush(stdout);
319 }
320
321 /* compute initial residual rt = r = F - J*dv */
322 J_times_dv (nvar, n1, n2, n3, dv, r, u);
323 for (int j = 0; j < ntotal; j++)
324 rt[j] = r[j] = F[j] - r[j];
325
326 *normres = norm2 (r, ntotal);
327 if (output == 1) {
328 printf ("bicgstab: %5d %10.3e\n", 0, (double) *normres);
329 fflush(stdout);
330 }
331
332 if (*normres <= tol)
333 return 0;
334
335 /* cgs iteration */
336 for (ii = 0; ii < itmax; ii++)
337 {
338 rho = scalarproduct (rt, r, ntotal);
339 if (fabs (rho) < rhotol)
340 break;
341
342 /* compute direction vector p */
343 if (ii == 0)
344 {
345
346 for (int j = 0; j < ntotal; j++)
347 p[j] = r[j];
348 }
349 else
350 {
351 beta = (rho / rho1) * (alpha / omega);
352
353 for (int j = 0; j < ntotal; j++)
354 p[j] = r[j] + beta * (p[j] - omega * vv[j]);
355 }
356
357 /* compute direction adjusting vector ph and scalar alpha */
358
359 for (int j = 0; j < ntotal; j++)
360 ph.d0[j] = 0;
361 for (int j = 0; j < NRELAX; j++) /* solves JFD*ph = p by relaxation*/
362 relax (ph.d0, nvar, n1, n2, n3, p, ncols, cols, JFD);
363
364 J_times_dv (nvar, n1, n2, n3, ph, vv, u); /* vv=J*ph*/
365 alpha = rho / scalarproduct (rt, vv, ntotal);
366
367 for (int j = 0; j < ntotal; j++)
368 s[j] = r[j] - alpha * vv[j];
369
370 /* early check of tolerance */
371 *normres = norm2 (s, ntotal);
372 if (*normres <= tol)
373 {
374
375 for (int j = 0; j < ntotal; j++)
376 dv.d0[j] += alpha * ph.d0[j];
377 if (output == 1) {
378 printf ("bicgstab: %5d %10.3e %10.3e %10.3e %10.3e\n",
379 ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
380 fflush(stdout);
381 }
382 break;
383 }
384
385 /* compute stabilizer vector sh and scalar omega */
386
387 for (int j = 0; j < ntotal; j++)
388 sh.d0[j] = 0;
389 for (int j = 0; j < NRELAX; j++) /* solves JFD*sh = s by relaxation*/
390 relax (sh.d0, nvar, n1, n2, n3, s, ncols, cols, JFD);
391
392 J_times_dv (nvar, n1, n2, n3, sh, t, u); /* t=J*sh*/
393 omega = scalarproduct (t, s, ntotal) / scalarproduct (t, t, ntotal);
394
395 /* compute new solution approximation */
396
397 for (int j = 0; j < ntotal; j++)
398 {
399 dv.d0[j] += alpha * ph.d0[j] + omega * sh.d0[j];
400 r[j] = s[j] - omega * t[j];
401 }
402 /* are we done? */
403 *normres = norm2 (r, ntotal);
404 if (output == 1) {
405 printf ("bicgstab: %5d %10.3e %10.3e %10.3e %10.3e\n",
406 ii + 1, (double) *normres, (double)alpha, (double)beta, (double)omega);
407 fflush(stdout);
408 }
409 if (*normres <= tol)
410 break;
411 rho1 = rho;
412 if (fabs (omega) < omegatol)
413 break;
414
415 }
416
417 /* free temporary storage */
418 free_dvector (r, 0, ntotal - 1);
419 free_dvector (p, 0, ntotal - 1);
420/* free_dvector(ph, 0, ntotal-1);*/
421 free_derivs (&ph, ntotal);
422 free_dvector (rt, 0, ntotal - 1);
423 free_dvector (s, 0, ntotal - 1);
424/* free_dvector(sh, 0, ntotal-1);*/
425 free_derivs (&sh, ntotal);
426 free_dvector (t, 0, ntotal - 1);
427 free_dvector (vv, 0, ntotal - 1);
428
429 free_dvector (F, 0, ntotal - 1);
430 free_derivs (&u, ntotal);
431
432 free_dmatrix (JFD, 0, ntotal - 1, 0, maxcol - 1);
433 free_imatrix (cols, 0, ntotal - 1, 0, maxcol - 1);
434 free_ivector (ncols, 0, ntotal - 1);
435
436 /* iteration failed */
437 if (ii > itmax)
438 return -1;
439
440 /* breakdown */
441 if (fabs (rho) < rhotol)
442 return -10;
443 if (fabs (omega) < omegatol)
444 return -11;
445
446 /* success! */
447 return ii + 1;
448}
449
450/* -------------------------------------------------------------------*/
451void
452TwoPunctures::Newton (int const nvar, int const n1, int const n2, int const n3,
453 derivs v,
454 double const tol, int const itmax)
455{
456
457
458 int ntotal = n1 * n2 * n3 * nvar, ii, it;
459 double *F, dmax, normres;
460 derivs u, dv;
461
462 F = dvector (0, ntotal - 1);
463 allocate_derivs (&dv, ntotal);
464 allocate_derivs (&u, ntotal);
465
466/* TestRelax(nvar, n1, n2, n3, v, dv.d0); */
467 it = 0;
468 dmax = 1;
469 while (dmax > tol && it < itmax)
470 {
471 if (it == 0)
472 {
473 F_of_v (nvar, n1, n2, n3, v, F, u);
474 dmax = norm_inf (F, ntotal);
475 }
476
477 for (int j = 0; j < ntotal; j++)
478 dv.d0[j] = 0;
479
480 if(verbose==1){
481 printf ("Newton: it=%d \t |F|=%e\n", it, (double)dmax);
482 printf ("bare mass: mp=%g \t mm=%g\n", (double) par_m_plus, (double) par_m_minus);
483 }
484
485 fflush(stdout);
486 ii =
487 bicgstab (nvar, n1, n2, n3, v, dv, verbose, 100, dmax * 1.e-3, &normres);
488
489 for (int j = 0; j < ntotal; j++)
490 v.d0[j] -= dv.d0[j];
491 F_of_v (nvar, n1, n2, n3, v, F, u);
492 dmax = norm_inf (F, ntotal);
493 it += 1;
494 }
495 if (itmax==0)
496 {
497 F_of_v (nvar, n1, n2, n3, v, F, u);
498 dmax = norm_inf (F, ntotal);
499 }
500
501 if(verbose==1)
502 printf ("Newton: it=%d \t |F|=%e \n", it, (double)dmax);
503
504 fflush(stdout);
505
506 free_dvector (F, 0, ntotal - 1);
507 free_derivs (&dv, ntotal);
508 free_derivs (&u, ntotal);
509}
510
511/* -------------------------------------------------------------------*/
512
513} // namespace TP
TP::TwoPunctures::free_derivs
void free_derivs(derivs *v, int n)
Definition FuncAndJacobian.cpp:69
TP::TwoPunctures::bicgstab
int bicgstab(int const nvar, int const n1, int const n2, int const n3, derivs v, derivs dv, int const output, int const itmax, double const tol, double *TP_RESTRICT const normres)
Definition Newton.cpp:277
TP::TwoPunctures::F_of_v
void F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F, derivs u)
Definition FuncAndJacobian.cpp:194
TP::TwoPunctures::Newton
void Newton(int nvar, int n1, int n2, int n3, derivs v, double tol, int itmax)
Definition Newton.cpp:452
TP::TwoPunctures::relax
void relax(double *TP_RESTRICT const dv, int const nvar, int const n1, int const n2, int const n3, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition Newton.cpp:176
TP::TwoPunctures::LineRelax_be
void LineRelax_be(double *TP_RESTRICT const dv, int const i, int const k, int const nvar, int const n1, int const n2, int const n3, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition Newton.cpp:117
TP::TwoPunctures::resid
void resid(double *TP_RESTRICT const res, int const ntotal, double const *TP_RESTRICT const dv, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition Newton.cpp:38
TP::TwoPunctures::LineRelax_al
void LineRelax_al(double *TP_RESTRICT const dv, int const j, int const k, int const nvar, int const n1, int const n2, int const n3, double const *TP_RESTRICT const rhs, int const *TP_RESTRICT const ncols, int const *TP_RESTRICT const *TP_RESTRICT const cols, double const *TP_RESTRICT const *TP_RESTRICT const JFD)
Definition Newton.cpp:57
TP::TwoPunctures::J_times_dv
void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u)
Definition FuncAndJacobian.cpp:372
TP::TwoPunctures::SetMatrix_JFD
void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix)
Definition FuncAndJacobian.cpp:594
TP::TwoPunctures::allocate_derivs
void allocate_derivs(derivs *v, int n)
Definition FuncAndJacobian.cpp:52
TP::TwoPunctures::norm_inf
double norm_inf(double const *TP_RESTRICT const F, int const ntotal)
Definition Newton.cpp:22
TP::TwoPunctures::Index
int Index(int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
Definition FuncAndJacobian.cpp:28
TP::TwoPunctures::TestRelax
void TestRelax(int nvar, int n1, int n2, int n3, derivs v, double *dv)
Definition Newton.cpp:225
TP::Utilities::norm2
double norm2(double *v, int n)
Definition TP_Utilities.cpp:506
TP::Utilities::dvector
double * dvector(long nl, long nh)
Definition TP_Utilities.cpp:30
TP::Utilities::free_dvector
void free_dvector(double *v, long nl, long nh)
Definition TP_Utilities.cpp:157
TP::Utilities::norm1
double norm1(double *v, int n)
Definition TP_Utilities.cpp:492
TP::Utilities::scalarproduct
double scalarproduct(double *v, double *w, int n)
Definition TP_Utilities.cpp:519
TP::Utilities::free_ivector
void free_ivector(int *v, long nl, long nh)
Definition TP_Utilities.cpp:149
TP::Utilities::ivector
int * ivector(long nl, long nh)
Definition TP_Utilities.cpp:16
TP::Utilities::dmatrix
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
Definition TP_Utilities.cpp:73
TP::Utilities::free_dmatrix
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
Definition TP_Utilities.cpp:174
TP::Utilities::free_imatrix
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
Definition TP_Utilities.cpp:165
TP::Utilities::imatrix
int ** imatrix(long nrl, long nrh, long ncl, long nch)
Definition TP_Utilities.cpp:44
TP
This file contains aliases for making access to the long state vector Q as used eg.
Definition CoordTransf.cpp:13
TP::derivs
struct TP::DERIVS derivs
TP::Step_Relax
const int Step_Relax
Definition Newton.cpp:18
TP::StencilSize
const int StencilSize
Definition Newton.cpp:15
TP::N_PlaneRelax
const int N_PlaneRelax
Definition Newton.cpp:16
TP::NRELAX
const int NRELAX
Definition Newton.cpp:17
TP::DERIVS::d0
double * d0
Definition TwoPunctures.h:31