File:	/home/sbrandt/cactus/Cactus/arrangements/CactusElliptic/EllSOR/src/Flat.c
1	: /*@@
2	:   @file      Flat.c
3	:   @date      Tue Aug 24 12:50:07 1999
4	:   @author    Gerd Lanfermann
5	:   @desc
6	:   SOR solver for 3D flat equation
7	:   @enddesc
8	: @@*/
9	:
10	:/*#define DEBUG_ELLIPTIC*/
11	:
12	:#include <stdlib.h>
13	:#include <stdio.h>
14	:#include <math.h>
15	:#include <string.h>
16	:
17	:#include "cctk.h"
18	:#include "cctk_Parameters.h"
19	:
20	:#include "Boundary.h"
21	:#include "Symmetry.h"
22	:#include "Ell_DBstructure.h"
23	:#include "EllBase.h"
24	:
25	:static const char *rcsid = "$Header$";
26	:
27	:CCTK_FILEVERSION(CactusElliptic_EllSOR_Flat_c)
28	:
29	:int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex,
30	:              CCTK_REAL *AbsTol, CCTK_REAL *RelTol);
31	:
32	:int SORFlat3D(cGH *GH, int FieldIndex, int MIndex, int NIndex,
33	:              CCTK_REAL *AbsTol, CCTK_REAL *RelTol)
34	:{
35	:  DECLARE_CCTK_PARAMETERS
36	:
37	:  int retval = ELL_SUCCESS;
38	:
39	:  /* The pointer to the data fields */
40	:  CCTK_REAL *Mlin=NULL;
41	:  CCTK_REAL *Nlin=NULL;
42	:  CCTK_REAL *var =NULL;
43	:
44	:  /* shortcuts for deltas,etc. */
45	:  CCTK_REAL dx,dy,dz;
46	:
47	:  /* Some physical variables */
48	:  CCTK_REAL omega, resnorm, residual, glob_residual, rjacobian;
49	:  CCTK_REAL finf;
50	:  CCTK_INT npow;
51	:  CCTK_REAL tol;
52	:
53	:  /* Iteration and stepping  variables */
54	:  int sorit;
55	:  CCTK_INT maxit;
56	:  int i, j, k;
57	:  int origin_sign;
58	:
59	:  /* stencil index */
60	:  int ijk;
61	:  int ipjk, ijpk, ijkp, imjk, ijmk, ijkm;
62	:
63	:  /* Coefficients for the stencil...  */
64	:  CCTK_REAL ac,ac_orig,aw,ae,an,as,at,ab;
65	:
66	:  /* Miscellaneous */
67	:  int sum_handle;
68	:  int sw[3];
69	:  int Mstorage=0, Nstorage=0;
70	:  static int firstcall = 1;
71	:  CCTK_REAL  dx2rec, dy2rec, dz2rec;
72	:  char *robin = NULL;
73	:  char *sor_accel = NULL;
74	:  int ierr = 0;
75	:  const void* input_array[1];
76	:  void*  reduction_value[1];
77	:  int input_array_dim[1];
78	:  int input_array_type_codes[1];
79	:
80	:  /* Boundary conditions */
81	:  int use_robin = 0;
82	:
83	:  input_array[0]             = (CCTK_REAL *)&resnorm;
84	:  input_array_dim[0]         = 0;
85	:  reduction_value[0]         = (CCTK_REAL *)&glob_residual;
86	:  input_array_type_codes[0]  = CCTK_VARIABLE_REAL;
87	:
88	:  /* Avoid compiler warnings */
89	:  RelTol = RelTol;
90	:
91	:  /* Get the reduction handle */
92	:  sum_handle = CCTK_LocalArrayReductionHandle("sum");
93	:  if (sum_handle<0)
94	:  {
95	:    CCTK_WARN(1,"SORFlat3D: Cannot get reduction handle for sum operation");
96	:  }
97	:
98	:  /* IF Robin BCs are set, prepare for a boundary call:
99	:     setup stencil width and get Robin constants (set by the routine
100	:     which is calling the solver interface) */
101	:
102	:  if (Ell_GetStrKey(&robin,"EllLinFlat::Bnd::Robin") < 0)
103	:  {
104	:    CCTK_WARN(2,"SORFlat3D: Robin keys not set");
105	:  }
106	:  else
107	:  {
108	:    if (CCTK_Equals(robin,"yes"))
109	:    {
110	:      use_robin = 1;
111	:
112	:      sw[0]=1;
113	:      sw[1]=1;
114	:      sw[2]=1;
115	:
116	:      if (Ell_GetRealKey(&finf, "EllLinFlat::Bnd::Robin::inf") ==
117	:          ELLGET_NOTSET)
118	:      {
119	:        CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::inf not set, "
120	:                  "setting to 1");
121	:        finf = 1;
122	:      }
123	:      if (Ell_GetIntKey(&npow, "EllLinFlat::Bnd::Robin::falloff")
124	:          == ELLGET_NOTSET)
125	:      {
126	:        CCTK_WARN(1,"SORFlat3D: EllLinFlat::Bnd::Robin::falloff not set, "
127	:                  "setting to 1");
128	:        npow = 1;
129	:      }
130	:    }
131	:  }
132	:
133	:  /* Get the maximum number of iterations allowed. */
134	:  if (Ell_GetIntKey(&maxit, "Ell::SORmaxit") == ELLGET_NOTSET)
135	:  {
136	:    CCTK_WARN(1,"SORFlat3D: Ell::SORmaxit not set. "
137	:              "Maximum allowed iterations being set to 100.");
138	:    maxit = 100;
139	:  }
140	:
141	:  /* Only supports absolute tolerance */
142	:  tol   = AbsTol[0];
143	:
144	:  /* Get the type of acceleration to use. */
145	:  if (Ell_GetStrKey(&sor_accel, "Ell::SORaccel") == ELLGET_NOTSET)
146	:  {
147	:    const char tmpstr[6] = "const";
148	:    CCTK_WARN(3, "SORFlat3D: Ell::SORaccel not set. "
149	:              "Omega being set to a constant value of 1.8.");
150	:    sor_accel = strdup(tmpstr);
151	:  }
152	:
153	:  /* Things to do only once! */
154	:  /* TODO: Need to handle this differently, since it may be called
155	:           for different reasons by the same code. */
156	:  if (firstcall==1)
157	:  {
158	:    if (CCTK_Equals(elliptic_verbose, "yes"))
159	:    {
160	:      if (CCTK_Equals(sor_accel, "cheb"))
161	:      {
162	:        CCTK_INFO("SOR with Chebyshev acceleration");
163	:      }
164	:      else if (CCTK_Equals(sor_accel, "const"))
165	:      {
166	:        CCTK_INFO("SOR with hardcoded omega = 1.8");
167	:      }
168	:      else if (CCTK_Equals(sor_accel, "none"))
169	:      {
170	:        CCTK_INFO("SOR with unaccelerated relaxation (omega = 1)");
171	:      }
172	:      else
173	:      {
174	:         CCTK_WARN(0, "SORFlat3D: sor_accel not set");
175	:      }
176	:    }
177	:    firstcall = 0;
178	:  }
179	:
180	:  /* Get the data ptr of these GFs, They all have to be
181	:     on the same timelevel; if we have a negative index for M/N,
182	:     this GF is not set,  there for don't even look for it and flag it  */
183	:
184	:  var = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, FieldIndex);
185	:  if (var == NULL)
186	:  {
187	:    CCTK_WARN(0, "SORFlat3D: Unable to get pointer to GF variable");
188	:  }
189	:
190	:  if (MIndex>=0)
191	:  {
192	:    Mlin = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, MIndex);
193	:    if (Mlin)
194	:    {
195	:      Mstorage = 1;
196	:    }
197	:    else
198	:    {
199	:      CCTK_WARN(0, "SORFlat3D: Unable to get pointer to M");
200	:    }
201	:  }
202	:
203	:  if (NIndex>=0)
204	:  {
205	:    Nlin = (CCTK_REAL *)CCTK_VarDataPtrI(GH, 0, NIndex);
206	:    if (Nlin)
207	:    {
208	:      Nstorage = 1;
209	:    }
210	:    else
211	:    {
212	:      CCTK_WARN(0, "SORFlat3D: Unable to get pointer to N");
213	:    }
214	:  }
215	:
216	:  /* Shortcuts */
217	:  dx   = GH->cctk_delta_space[0];
218	:  dy   = GH->cctk_delta_space[1];
219	:  dz   = GH->cctk_delta_space[2];
220	:
221	:  dx2rec = 1.0/(dx*dx);
222	:  dy2rec = 1.0/(dy*dy);
223	:  dz2rec = 1.0/(dz*dz);
224	:
225	:  ae = dx2rec;
226	:  aw = dx2rec;
227	:  an = dy2rec;
228	:  as = dy2rec;
229	:  at = dz2rec;
230	:  ab = dz2rec;
231	:
232	:  ac_orig = -2.0*dx2rec - 2.0*dy2rec - 2.0*dz2rec;
233	:
234	:  /* Initialize omega. */
235	:  /* TODO: Make it so that the value of the constant omega can be set. */
236	:  if (CCTK_Equals(sor_accel, "const"))
237	:  {
238	:    omega = 1.8;
239	:  }
240	:  else
241	:  {
242	:    omega = 1.;
243	:  }
244	:
245	:  /* Set the spectral radius of the Jacobi iteration. */
246	:  /* TODO: I think this can be easily computed for flat metrics? */
247	:  rjacobian =  0.999;
248	:
249	:  /* Compute whether the origin of this processor's sub grid is "red" or
250	:     "black". */
251	:  origin_sign = (GH->cctk_lbnd[0] + GH->cctk_lbnd[1] + GH->cctk_lbnd[2])%2;
252	:
253	:  /* start at 1 for historic (Fortran) reasons */
254	:  for (sorit=1; sorit<=maxit; sorit++)
255	:  {
256	:    resnorm = 0.;
257	:
258	:    for (k=1; k<GH->cctk_lsh[2]-1; k++)
259	:    {
260	:      for (j=1; j<GH->cctk_lsh[1]-1; j++)
261	:      {
262	:        for (i=1; i<GH->cctk_lsh[0]-1; i++)
263	:        {
264	:          if ((origin_sign + i + j + k)%2 == sorit%2)
265	:          {
266	:            ac = ac_orig;
267	:
268	:            ijk   = CCTK_GFINDEX3D(GH,i  ,j  ,k  );
269	:            ipjk  = CCTK_GFINDEX3D(GH,i+1,j  ,k  );
270	:            imjk  = CCTK_GFINDEX3D(GH,i-1,j  ,k  );
271	:            ijpk  = CCTK_GFINDEX3D(GH,i  ,j+1,k  );
272	:            ijmk  = CCTK_GFINDEX3D(GH,i  ,j-1,k  );
273	:            ijkp  = CCTK_GFINDEX3D(GH,i  ,j  ,k+1);
274	:            ijkm  = CCTK_GFINDEX3D(GH,i  ,j  ,k-1);
275	:
276	:            if (Mstorage)
277	:            {
278	:              ac += Mlin[ijk];
279	:            }
280	:
281	:            residual = ac * var[ijk]
282	:                + ae*var[ipjk] + aw*var[imjk]
283	:                + an*var[ijpk] + as*var[ijmk]
284	:                + at*var[ijkp] + ab*var[ijkm];
285	:
286	:            if (Nstorage)
287	:            {
288	:              residual += Nlin[ijk];
289	:            }
290	:
291	:            resnorm += fabs(residual);
292	:
293	:            var[ijk] -= omega*residual/ac;
294	:          }
295	:        }
296	:      }
297	:    }
298	:
299	:    /* reduction operation on processor-local residual values */
300	:    ierr = CCTK_ReduceArraysGlobally(GH, -1,sum_handle, -1, 1,input_array,0,
301	:                          input_array_dim,
302	:                          input_array_type_codes,
303	:                          1,
304	:                          input_array_type_codes,
305	:                          reduction_value);
306	:
307	:    if (ierr<0)
308	:    {
309	:      CCTK_WARN(1,"SORFlat3D: Reduction of Norm failed");
310	:    }
311	:
312	:#ifdef DEBUG_ELLIPTIC
313	:    printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol);
314	:#endif
315	:
316	:    glob_residual /= (GH->cctk_gsh[0]*GH->cctk_gsh[1]*GH->cctk_gsh[2]);
317	:
318	:
319	:#ifdef DEBUG_ELLIPTIC
320	:    printf("it: %d SOR solve residual %f (tol %f)\n",sorit,glob_residual,tol);
321	:#endif
322	:
323	:    /* apply symmetry boundary conditions within loop */
324	:    if (CartSymVI(GH,FieldIndex)<0)
325	:    {
326	:      CCTK_WARN(1,"SORFlat3D: CartSymVI failed in EllSOR loop");
327	:      break;
328	:    }
329	:
330	:    /* apply Robin boundary conditions within loop */
331	:    if (use_robin)
332	:    {
333	:      if (BndRobinVI(GH, sw, finf, npow,  FieldIndex)<0)
334	:      {
335	:        CCTK_WARN(1,"SORFlat3D: BndRobinVI failed in EllSOR loop");
336	:        break;
337	:      }
338	:    }
339	:
340	:    /* synchronization of grid variable */
341	:    CCTK_SyncGroupWithVarI(GH, FieldIndex);
342	:
343	:    /* Leave iteration loop if tolerance criterium is met */
344	:    if (glob_residual<tol)
345	:    {
346	:      break;
347	:    }
348	:
349	:    /* Update omega if using Chebyshev acceleration. */
350	:    if (CCTK_Equals(sor_accel, "cheb"))
351	:    {
352	:      if (sorit == 1)
353	:      {
354	:        omega = 1. / (1. - 0.5 * rjacobian * rjacobian);
355	:      }
356	:      else
357	:      {
358	:         omega = 1. / (1. - 0.25 * rjacobian * rjacobian * omega);
359	:      }
360	:    }
361	:
362	:  }
363	:
364	:  if (elliptic_verbose)
365	:  {
366	:    CCTK_VInfo("EllSOR","Required SOR tolerence was set at %f",tol);
367	:    CCTK_VInfo("EllSOR","Achieved SOR residual was %f",glob_residual);
368	:    CCTK_VInfo("EllSOR","Number of iterations was %d (max: %d)",sorit,maxit);
369	:  }
370	:
371	:  if (glob_residual>tol)
372	:  {
373	:    CCTK_WARN(1,"SORFlat3D: SOR Solver did not converge");
374	:    retval = ELL_NOCONVERGENCE;
375	:  }
376	:
377	:  if (use_robin)
378	:  {
379	:    free(robin);
380	:  }
381	:  if (sor_accel)
382	:  {
383	:    free(sor_accel);
384	:  }
385	:
386	:  return retval;
387	:}