Peano 4
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>
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
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
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
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
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ our matrix elements will nabla phi_i dx f By this will be a sparse as these basis functions are chosen to not overlap with each other almost everywhere In other they have only local support We can read off the right hand side by taking our known right hand side f$ f f$ and integrating against an appropriate test phi_i dx f Please excuse the slight abuse of notation here There should probably be a clearer indication that we move from a continuous f$ f f$ to some discrete f$ f_i f$ We can demonstrate simply It s worth as when we discuss the discontinuous version of this it will no longer disappear We take our left hand side and discretise it
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ j
we integrate over each cell and then take the sum across each of the cells We also consider the terms that enter the f$ k
#define TP_RESTRICT
We compile TwoPunctures-Standalone as C++, and C++ does not yet now about the C99 language feature "r...
void free_derivs(derivs *v, int n)
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
void F_of_v(int nvar, int n1, int n2, int n3, derivs v, double *F, derivs u)
void Newton(int nvar, int n1, int n2, int n3, derivs v, double tol, int itmax)
Definition Newton.cpp:452
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
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
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
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
void J_times_dv(int nvar, int n1, int n2, int n3, derivs dv, double *Jdv, derivs u)
void SetMatrix_JFD(int nvar, int n1, int n2, int n3, derivs u, int *ncols, int **cols, double **Matrix)
void allocate_derivs(derivs *v, int n)
double norm_inf(double const *TP_RESTRICT const F, int const ntotal)
Definition Newton.cpp:22
int Index(int ivar, int i, int j, int k, int nvar, int n1, int n2, int n3)
void TestRelax(int nvar, int n1, int n2, int n3, derivs v, double *dv)
Definition Newton.cpp:225
examples::exahype2::elastic::VariableShortcuts s
Definition loh.cpp:10
double f(const tarch::la::Vector< Dimensions, double > &x)
double norm2(double *v, int n)
double * dvector(long nl, long nh)
void free_dvector(double *v, long nl, long nh)
double norm1(double *v, int n)
double scalarproduct(double *v, double *w, int n)
void free_ivector(int *v, long nl, long nh)
int * ivector(long nl, long nh)
double ** dmatrix(long nrl, long nrh, long ncl, long nch)
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
int ** imatrix(long nrl, long nrh, long ncl, long nch)
This file contains aliases for making access to the long state vector Q as used eg.
const int Step_Relax
Definition Newton.cpp:18
const int StencilSize
Definition Newton.cpp:15
const int N_PlaneRelax
Definition Newton.cpp:16
const int NRELAX
Definition Newton.cpp:17
examples::exahype2::elastic::VariableShortcuts sh
double * d0