1 /* PoissonMG.c */
2 /**********************************************************************************************************
3 Copyright (c) 2002-2013 Abdul-Rahman Allouche. All rights reserved
4
5 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated
6 documentation files (the Gabedit), to deal in the Software without restriction, including without limitation
7 the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software,
8 and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
9
10 The above copyright notice and this permission notice shall be included in all copies or substantial portions
11 of the Software.
12
13 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
14 TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
15 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
16 CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
17 DEALINGS IN THE SOFTWARE.
18 ************************************************************************************************************/
19
20 #include "../../Config.h"
21 #include <math.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include <stdlib.h>
25 #include <ctype.h>
26 #include <gtk/gtk.h>
27 #include <glib/gi18n.h>
28 #include "../Utils/Vector3d.h"
29 #include "../Utils/Transformation.h"
30 #include "../Utils/Constants.h"
31 #include "../Utils/Zlm.h"
32 #include "../Utils/MathFunctions.h"
33 #include "PoissonMG.h"
34 #include "../Common/GabeditType.h"
35 #include "../Display/GlobalOrb.h"
36 #include "../Display/StatusOrb.h"
37
38 /*********************************************************/
getPoissonUsingDomain(DomainMG * domain)39 PoissonMG* getPoissonUsingDomain(DomainMG* domain)
40 {
41 PoissonMG * ps = g_malloc(sizeof(PoissonMG));
42
43 ps->potential = getNewGridMGUsingDomain(domain);
44 ps->source = getNewGridMGUsingDomain(domain);
45
46 setOperationGridMG(ps->potential, GABEDIT_ALL);
47 setOperationGridMG(ps->source, GABEDIT_ALL);
48
49 ps->condition = GABEDIT_CONDITION_MULTIPOL;
50 ps->setBoundary = NULL;
51
52 ps->diag = 1;
53 return ps;
54 }
55 /*********************************************************/
getNullPoissonMG()56 PoissonMG* getNullPoissonMG()
57 {
58 DomainMG domain = getNullDomainMG();
59 return getPoissonUsingDomain(&domain);
60 }
61 /*********************************************************/
getPoissonMG2(GridMG * p,GridMG * s,Condition c,SetBoundary set)62 PoissonMG* getPoissonMG2(GridMG* p, GridMG* s, Condition c, SetBoundary set)
63 {
64 if(!s) return getNullPoissonMG();
65 /*
66 PoissonMG* ps = getPoissonUsingDomain(&(s->domain));
67 if(!p) p = getNewGridMGUsingDomain(&(s->domain));
68 */
69 PoissonMG * ps = g_malloc(sizeof(PoissonMG));
70 ps->potential = p;
71 ps->source = s;
72 ps->diag = getDiagGridMG(ps->potential);
73
74 ps->condition = c;
75 ps->setBoundary = NULL;
76 return ps;
77 }
78 /*********************************************************/
getPoissonMG(GridMG * p,GridMG * s)79 PoissonMG* getPoissonMG(GridMG* p, GridMG* s)
80 {
81 return getPoissonMG2(p,s, GABEDIT_CONDITION_MULTIPOL, NULL);
82 }
83 /*********************************************************/
getCopyPoissonMG(PoissonMG * ps)84 PoissonMG* getCopyPoissonMG(PoissonMG* ps)
85 {
86 if(!ps) return getNullPoissonMG();
87
88 PoissonMG * newps = g_malloc(sizeof(PoissonMG));
89 newps->potential = getNewGridMGFromOldGrid(ps->potential);
90 newps->source = getNewGridMGFromOldGrid(ps->source);
91 newps->diag = ps->diag ;
92
93 newps->condition = ps->condition ;
94 newps->setBoundary = ps->setBoundary;
95 return newps;
96 }
97 /*********************************************************/
destroyPoissonMG(PoissonMG * ps)98 void destroyPoissonMG(PoissonMG* ps)
99 {
100 destroyGridMG(ps->potential);
101 destroyGridMG(ps->source);
102 ps->diag = 1;
103 }
104 /*********************************************************/
setOperationPoissonMG(PoissonMG * ps,OperationTypeMG operation)105 void setOperationPoissonMG(PoissonMG* ps, OperationTypeMG operation)
106 {
107 setOperationGridMG(ps->potential, operation);
108 setOperationGridMG(ps->source, operation);
109 }
110 /*********************************************************/
getDomainPoissonMG(PoissonMG * ps)111 DomainMG getDomainPoissonMG(PoissonMG* ps)
112 {
113 return getDomainGridMG(ps->potential);
114 }
115 /*********************************************************/
getDiagPoissonMG(PoissonMG * ps)116 gdouble getDiagPoissonMG(PoissonMG* ps)
117 {
118 return getDiagGridMG(ps->potential);
119 }
120 /*********************************************************/
prolongationPoissonMG(PoissonMG * ps)121 void prolongationPoissonMG(PoissonMG* ps)
122 {
123 prolongationGridMG(ps->potential);
124 levelUpGridMG(ps->source);
125 ps->diag = getDiagGridMG(ps->potential);
126 }
127 /*********************************************************/
interpolationTriLinearPoissonMG(PoissonMG * ps)128 void interpolationTriLinearPoissonMG(PoissonMG* ps)
129 {
130 interpolationTriLinearGridMG(ps->potential);
131 levelUpGridMG(ps->source);
132 ps->diag = getDiagGridMG(ps->potential);
133 }
134 /*********************************************************/
interpolationCubicPoissonMG(PoissonMG * ps)135 void interpolationCubicPoissonMG(PoissonMG* ps)
136 {
137 interpolationCubicGridMG(ps->potential);
138 levelUpGridMG(ps->source);
139 ps->diag = getDiagGridMG(ps->potential);
140 }
141 /*********************************************************/
restrictionPoissonMG(PoissonMG * ps)142 void restrictionPoissonMG(PoissonMG* ps)
143 {
144 DomainMG domain = getDomainPoissonMG(ps);
145 GridMG* tau = getNewGridMGUsingDomain(&domain);
146 laplacianGridMG(tau,ps->potential);
147
148 restrictionGridMG(ps->potential);
149 /*restrictionInjectionGridMG(ps->potential);*/
150 tradesBoundaryPoissonMG(ps); /* Important pour le calcul du laplacien*/
151 restrictionGridMG(ps->source);
152 restrictionGridMG(tau);
153
154 multEqualRealGridMG(tau,-1.0);
155 ps->diag = plusLaplacianGridMG(tau,ps->potential);
156
157 plusEqualGridMG(ps->source,tau);
158 }
159 /*********************************************************/
restrictionInjectionPoissonMG(PoissonMG * ps)160 void restrictionInjectionPoissonMG(PoissonMG* ps)
161 {
162 DomainMG domain = getDomainPoissonMG(ps);
163 GridMG* tau = getNewGridMGUsingDomain(&domain);
164 laplacianGridMG(tau,ps->potential);
165
166 restrictionInjectionGridMG(ps->potential);
167 tradesBoundaryPoissonMG(ps); /* Important pour le calcul du laplacien*/
168 restrictionInjectionGridMG(ps->source);
169 restrictionInjectionGridMG(tau);
170
171 multEqualRealGridMG(tau,-1.0);
172 ps->diag = plusLaplacianGridMG(tau,ps->potential);
173
174 plusEqualGridMG(ps->source,tau);
175 destroyGridMG(tau);
176 }
177 /*********************************************************/
residualPoissonMG(PoissonMG * ps)178 GridMG* residualPoissonMG(PoissonMG* ps)
179 {
180 OperationTypeMG opPotential = getOperationGridMG(ps->potential);
181 OperationTypeMG opSource = getOperationGridMG(ps->source);
182 GridMG* res=getNewGridMG();
183
184 setOperationGridMG(res, GABEDIT_INTERIOR);
185 setOperationGridMG(ps->potential, GABEDIT_INTERIOR);
186 setOperationGridMG(ps->source, GABEDIT_INTERIOR);
187
188 copyGridMG(res,ps->source);
189 moinsLaplacianGridMG(res,ps->potential);
190
191 setOperationGridMG(res, GABEDIT_BOUNDARY);
192 initGridMG(res,0.0);
193
194 setOperationGridMG(ps->potential, opPotential);
195 setOperationGridMG(ps->source, opSource);
196 setOperationGridMG(res, GABEDIT_INTERIOR);
197
198 return res;
199 }
200 /*********************************************************/
setBoundaryEwaldPoissonMG(PoissonMG * ps)201 void setBoundaryEwaldPoissonMG(PoissonMG* ps)
202 {
203 gdouble x,y,z;
204 gdouble xs,ys,zs;
205 gdouble Q = 0;
206 int ixs,iys,izs;
207 gdouble xOff =0, yOff = 0, zOff = 0;
208 DomainMG domain = getDomainPoissonMG(ps);
209 gdouble v = -domain.cellVolume/4/PI;
210 static gdouble PRECISION = 1e-12;
211
212 printf(_("Set boundaries using ewald sum\n"));
213 setOperationGridMG(ps->potential,GABEDIT_BOUNDARY);
214 initGridMG(ps->potential , 0.0);
215 setOperationGridMG(ps->potential,GABEDIT_ALL);
216
217 for(ixs=domain.iXBeginInterior;ixs<=domain.iXEndInterior;ixs++)
218 for(iys = domain.iYBeginInterior;iys <=domain.iYEndInterior;iys++)
219 for(izs = domain.iZBeginInterior;izs <=domain.iZEndInterior;izs++)
220 {
221 xs = domain.x0 + ixs*domain.xh - xOff;
222 ys = domain.y0 + iys*domain.yh - yOff;
223 zs = domain.z0 + izs*domain.zh - zOff;
224 Q = getValGridMG(ps->source, ixs, iys, izs)*v;
225
226 int ix,iy,iz;
227
228 gdouble x0 = domain.x0 - xOff-xs;
229 gdouble y0 = domain.y0 - yOff-ys;
230 gdouble z0 = domain.z0 - zOff-zs;
231 gdouble invR2;
232
233 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
234 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
235 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
236 {
237 x = x0 + ix*domain.xh ;
238 y = y0 + iy*domain.yh ;
239 z = z0 + iz*domain.zh ;
240
241 invR2 = 1.0 / (x*x + y*y + z*z + PRECISION);
242 addValGridMG(ps->potential, ix, iy, iz, sqrt(invR2) *Q);
243 }
244 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
245 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
246 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
247 {
248 x = x0 + ix*domain.xh ;
249 y = y0 + iy*domain.yh ;
250 z = z0 + iz*domain.zh ;
251
252 invR2 = 1.0 / (x*x + y*y + z*z + PRECISION);
253 addValGridMG(ps->potential, ix, iy, iz, sqrt(invR2) *Q);
254 }
255
256 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
257 {
258 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
259 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
260 {
261 x = x0 + ix*domain.xh ;
262 y = y0 + iy*domain.yh ;
263 z = z0 + iz*domain.zh ;
264
265 invR2 = 1.0 / (x*x + y*y + z*z + PRECISION);
266 addValGridMG(ps->potential, ix, iy, iz, sqrt(invR2) *Q);
267 }
268 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
269 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
270 {
271 x = x0 + ix*domain.xh ;
272 y = y0 + iy*domain.yh ;
273 z = z0 + iz*domain.zh ;
274
275 invR2 = 1.0 / (x*x + y*y + z*z + PRECISION);
276 addValGridMG(ps->potential, ix, iy, iz, sqrt(invR2) *Q);
277 }
278 }
279 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
280 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
281 {
282 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
283 {
284 x = x0 + ix*domain.xh ;
285 y = y0 + iy*domain.yh ;
286 z = z0 + iz*domain.zh ;
287
288 invR2 = 1.0 / (x*x + y*y + z*z + PRECISION);
289 addValGridMG(ps->potential, ix, iy, iz, sqrt(invR2) *Q);
290 }
291 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
292 {
293 x = x0 + ix*domain.xh ;
294 y = y0 + iy*domain.yh ;
295 z = z0 + iz*domain.zh ;
296
297 invR2 = 1.0 / (x*x + y*y + z*z + PRECISION);
298 addValGridMG(ps->potential, ix, iy, iz, sqrt(invR2) *Q);
299 }
300 }
301 }
302 }
303 /*********************************************************/
getCOff(PoissonMG * ps,gdouble * pxOff,gdouble * pyOff,gdouble * pzOff)304 static void getCOff(PoissonMG* ps, gdouble* pxOff, gdouble* pyOff, gdouble* pzOff)
305 {
306 gdouble temp;
307 gdouble x,y,z;
308 DomainMG domain = getDomainPoissonMG(ps);
309 int ixs,iys,izs;
310 gdouble xOff =0, yOff = 0, zOff = 0;
311 gdouble Q = 0;
312
313 for(ixs=domain.iXBeginInterior;ixs<=domain.iXEndInterior;ixs++)
314 for(iys = domain.iYBeginInterior;iys <=domain.iYEndInterior;iys++)
315 for(izs = domain.iZBeginInterior;izs <=domain.iZEndInterior;izs++)
316 {
317 x = domain.x0 + ixs*domain.xh;
318 y = domain.y0 + iys*domain.yh;
319 z = domain.z0 + izs*domain.zh;
320 temp = getValGridMG(ps->source,ixs, iys, izs);
321 Q += temp;
322 xOff += temp*x;
323 yOff += temp*y;
324 zOff += temp*z;
325 }
326 if(Q!=0)
327 {
328 *pxOff = xOff/Q;
329 *pyOff = yOff/Q;
330 *pzOff = zOff/Q;
331 }
332 else
333 {
334 *pxOff = 0;
335 *pyOff = 0;
336 *pzOff = 0;
337 }
338 }
339 /*********************************************************/
setBoundaryMultipolPoissonMG(PoissonMG * ps)340 void setBoundaryMultipolPoissonMG(PoissonMG* ps)
341 {
342 const int lmax = 3;
343 Zlm zlm[lmax+1][2*lmax+1+1];
344 gint l,m;
345 gint i;
346 static gdouble PRECISION = 1e-12;
347 gdouble temp;
348 gdouble x,y,z;
349 DomainMG domain = getDomainPoissonMG(ps);
350 gdouble Q[lmax+1][2*lmax+1+1];
351 int ixs,iys,izs;
352 gdouble v = -domain.cellVolume/4/PI;
353 gdouble R[3];
354 gdouble xOff =0, yOff = 0, zOff = 0;
355 gdouble r;
356
357 /* printf("Begin Set boundaries using multipole approximation\n");*/
358
359 for( l=0; l<=lmax; l++)
360 for( m=-l; m<=l; m++)
361 {
362 zlm[l][m+l] = getZlm(l, m);
363 }
364 /* printf("End Zlm\n");*/
365
366 for( l=0; l<=lmax; l++)
367 for( m=-l; m<=l; m++)
368 Q[l][m+l] = 0.0;
369
370
371
372 printf(_("Set boundaries using multipole approximation\n"));
373
374 setOperationGridMG(ps->potential,GABEDIT_BOUNDARY);
375 initGridMG(ps->potential,0.0);
376 /* printf("End initGridMG\n");*/
377 setOperationGridMG(ps->potential,GABEDIT_ALL);
378 getCOff(ps, &xOff, &yOff, &zOff);
379
380 for(ixs=domain.iXBeginInterior;ixs<=domain.iXEndInterior;ixs++)
381 for(iys = domain.iYBeginInterior;iys <=domain.iYEndInterior;iys++)
382 for(izs = domain.iZBeginInterior;izs <=domain.iZEndInterior;izs++)
383 {
384 x = domain.x0 + ixs*domain.xh - xOff;
385 y = domain.y0 + iys*domain.yh - yOff;
386 z = domain.z0 + izs*domain.zh - zOff;
387 temp = getValGridMG(ps->source,ixs, iys, izs)*v;
388 r = sqrt(x*x + y*y + z*z+PRECISION);
389 R[0] = x;
390 R[1] = y;
391 R[2] = z;
392 if(r>0)
393 for( i=0; i<3; i++)
394 R[i] /= r;
395
396 for( l=0; l<=lmax; l++)
397 {
398 gdouble p = temp*pow(r,l);
399 for( m=-l; m<=l; m++)
400 {
401 Q[l][m+l] += p*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2]);
402 }
403 }
404 }
405 printf(_("Total charge = %f\n"),Q[0][0]*sqrt(4*PI));
406 printf(_("Center = %f %f %f\n"),xOff, yOff, zOff);
407 for( l=0; l<=lmax; l++)
408 for( m=-l; m<=l; m++)
409 {
410
411 unsigned int absm = abs(m);
412 gdouble Norm = 1;
413 Norm = sqrt((2*l+1)/(4*PI))*sqrt(factorial(l+absm)*factorial(l-absm))/factorial(l)/pow(2.0,absm);
414 if(m!=0) Norm *= sqrt(2.0);
415 printf("Q(%d,%d)=%f\n",l,m,Q[l][m+l]/Norm);
416 Q[l][m+l] *= 4*PI/(2*l+1);
417 }
418 int ix,iy,iz;
419
420 gdouble x0 = domain.x0 - xOff;
421 gdouble y0 = domain.y0 - yOff;
422 gdouble z0 = domain.z0 - zOff;
423 gdouble invR;
424
425 for(ix=domain.iXBeginBoundaryLeft;ix<=domain.iXEndBoundaryLeft;ix++)
426 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
427 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
428 {
429 x = x0 + ix*domain.xh ;
430 y = y0 + iy*domain.yh ;
431 z = z0 + iz*domain.zh ;
432
433 invR = 1.0 /sqrt (x*x + y*y + z*z + PRECISION);
434 R[0] = x*invR;
435 R[1] = y*invR;
436 R[2] = z*invR;
437 gdouble v = 0;
438 for( l=0; l<=lmax; l++)
439 {
440 temp = pow(invR,l+1);
441 for( m=-l; m<=l; m++)
442 {
443 if(fabs(Q[l][m+l])<10*PRECISION) continue;
444 v += temp*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2])*Q[l][m+l];
445 }
446 }
447 setValGridMG(ps->potential, ix, iy, iz, v);
448 }
449 for(ix=domain.iXBeginBoundaryRight;ix<=domain.iXEndBoundaryRight;ix++)
450 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
451 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
452 {
453 x = x0 + ix*domain.xh ;
454 y = y0 + iy*domain.yh ;
455 z = z0 + iz*domain.zh ;
456
457 invR = 1.0 /sqrt (x*x + y*y + z*z + PRECISION);
458 R[0] = x*invR;
459 R[1] = y*invR;
460 R[2] = z*invR;
461 gdouble v = 0;
462 for( l=0; l<=lmax; l++)
463 {
464 temp = pow(invR,l+1);
465 for( m=-l; m<=l; m++)
466 {
467 if(fabs(Q[l][m+l])<10*PRECISION) continue;
468 v += temp*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2])*Q[l][m+l];
469 }
470 }
471 setValGridMG(ps->potential, ix, iy, iz, v);
472 }
473
474 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
475 {
476 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryLeft;iy++)
477 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
478 {
479 x = x0 + ix*domain.xh ;
480 y = y0 + iy*domain.yh ;
481 z = z0 + iz*domain.zh ;
482 invR = 1.0 /sqrt (x*x + y*y + z*z + PRECISION);
483 R[0] = x*invR;
484 R[1] = y*invR;
485 R[2] = z*invR;
486 gdouble v = 0;
487 for( l=0; l<=lmax; l++)
488 {
489 temp = pow(invR,l+1);
490 for( m=-l; m<=l; m++)
491 {
492 if(fabs(Q[l][m+l])<10*PRECISION) continue;
493 v += temp*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2])*Q[l][m+l];
494 }
495 }
496 setValGridMG(ps->potential, ix, iy, iz, v);
497
498 }
499 for(iy = domain.iYBeginBoundaryRight;iy <=domain.iYEndBoundaryRight;iy++)
500 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryRight;iz++)
501 {
502 x = x0 + ix*domain.xh ;
503 y = y0 + iy*domain.yh ;
504 z = z0 + iz*domain.zh ;
505 invR = 1.0 /sqrt (x*x + y*y + z*z + PRECISION);
506 R[0] = x*invR;
507 R[1] = y*invR;
508 R[2] = z*invR;
509 gdouble v = 0;
510 for( l=0; l<=lmax; l++)
511 {
512 temp = pow(invR,l+1);
513 for( m=-l; m<=l; m++)
514 {
515 if(fabs(Q[l][m+l])<10*PRECISION) continue;
516 v += temp*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2])*Q[l][m+l];
517 }
518 }
519 setValGridMG(ps->potential, ix, iy, iz, v);
520
521 }
522 }
523 for(ix = domain.iXBeginBoundaryLeft;ix <=domain.iXEndBoundaryRight;ix++)
524 for(iy = domain.iYBeginBoundaryLeft;iy <=domain.iYEndBoundaryRight;iy++)
525 {
526 for(iz = domain.iZBeginBoundaryLeft;iz <=domain.iZEndBoundaryLeft;iz++)
527 {
528 x = x0 + ix*domain.xh ;
529 y = y0 + iy*domain.yh ;
530 z = z0 + iz*domain.zh ;
531
532 invR = 1.0 /sqrt (x*x + y*y + z*z + PRECISION);
533 R[0] = x*invR;
534 R[1] = y*invR;
535 R[2] = z*invR;
536 gdouble v = 0;
537 for( l=0; l<=lmax; l++)
538 {
539 temp = pow(invR,l+1);
540 for( m=-l; m<=l; m++)
541 {
542 if(fabs(Q[l][m+l])<10*PRECISION) continue;
543 v += temp*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2])*Q[l][m+l];
544 }
545 }
546 setValGridMG(ps->potential, ix, iy, iz, v);
547 }
548 for(iz = domain.iZBeginBoundaryRight;iz <=domain.iZEndBoundaryRight;iz++)
549 {
550 x = x0 + ix*domain.xh ;
551 y = y0 + iy*domain.yh ;
552 z = z0 + iz*domain.zh ;
553
554 invR = 1.0 /sqrt (x*x + y*y + z*z + PRECISION);
555 R[0] = x*invR;
556 R[1] = y*invR;
557 R[2] = z*invR;
558 gdouble v = 0;
559 for( l=0; l<=lmax; l++)
560 {
561 temp = pow(invR,l+1);
562 for( m=-l; m<=l; m++)
563 {
564 if(fabs(Q[l][m+l])<10*PRECISION) continue;
565 v += temp*getValueZlm(&zlm[l][m+l],R[0],R[1],R[2])*Q[l][m+l];
566 }
567 }
568 setValGridMG(ps->potential, ix, iy, iz, v);
569
570 }
571 }
572 }
573 /*********************************************************/
residualNormPoissonMG(PoissonMG * ps)574 gdouble residualNormPoissonMG(PoissonMG* ps)
575 {
576 GridMG* res;
577 gdouble normal = 1;
578
579 res = residualPoissonMG(ps);
580 setOperationGridMG(res, GABEDIT_INTERIOR);
581 normal = normGridMG(res);
582 destroyGridMG(res);
583
584 return normal;
585 }
586 /*********************************************************/
tradesBoundaryPoissonMG(PoissonMG * ps)587 void tradesBoundaryPoissonMG(PoissonMG* ps)
588 {
589 switch(ps->condition)
590 {
591 case GABEDIT_CONDITION_PERIODIC :
592 case GABEDIT_CONDITION_CLUSTER :
593 tradesBoundaryGridMG(ps->potential, ps->condition);
594 break;
595 case GABEDIT_CONDITION_MULTIPOL :
596 setBoundaryMultipolPoissonMG(ps);
597 break;
598 case GABEDIT_CONDITION_EWALD :
599 setBoundaryEwaldPoissonMG(ps);
600 break;
601 case GABEDIT_CONDITION_EXTERNAL :
602 ps->setBoundary(ps->potential, ps->source);
603 break;
604 }
605 }
606 /*********************************************************/
smootherPoissonMG(PoissonMG * ps,int max)607 void smootherPoissonMG(PoissonMG* ps, int max)
608 {
609 OperationTypeMG opPotential = getOperationGridMG(ps->potential);
610 OperationTypeMG opSource = getOperationGridMG(ps->source);
611 GridMG* res =getNewGridMG();
612 gint i;
613 gdouble rms;
614
615 setOperationGridMG(res, GABEDIT_INTERIOR);
616 setOperationGridMG(ps->potential, GABEDIT_INTERIOR);
617 setOperationGridMG(ps->source, GABEDIT_INTERIOR);
618
619 for(i=0;i<max;i++)
620 {
621 if(ps->condition == GABEDIT_CONDITION_PERIODIC) tradesBoundaryPoissonMG(ps);
622 res = residualPoissonMG(ps);
623 rms = normGridMG(res);
624 /* printf("rms = %f\n",rms);*/
625 multEqualRealGridMG(res, ps->diag*0.8);
626 plusEqualGridMG(ps->potential, res);
627 destroyGridMG(res);
628 }
629 setOperationGridMG(ps->potential, opPotential);
630 setOperationGridMG(ps->source, opSource);
631 }
632 /*********************************************************/
printFilePoissonMG(PoissonMG * ps)633 void printFilePoissonMG(PoissonMG* ps)
634 {
635 FILE* file;
636 file = fopen("Poisson.cube","w");
637 printGridFileGridMG(ps->potential, file);
638 fclose(file);
639 }
640 /*********************************************************/
printFileNamePoissonMG(PoissonMG * ps,char * fileName)641 void printFileNamePoissonMG(PoissonMG* ps, char* fileName)
642 {
643 FILE* file;
644
645 file = fopen(fileName,"w");
646 if(!file) return;
647
648 printGridFileGridMG(ps->potential, file);
649 fclose(file);
650 }
651 /*********************************************************/
printMinPoissonMG(PoissonMG * ps)652 void printMinPoissonMG(PoissonMG* ps)
653 {
654 printf("Potential\t");
655 printMinGridMG(ps->potential);
656 printf("Source\t");
657 printMinGridMG(ps->source);
658 }
659 /*********************************************************/
printMaxPoissonMGAll(PoissonMG * ps)660 void printMaxPoissonMGAll(PoissonMG* ps)
661 {
662 printf("Potential\t");
663 printMaxGridMG(ps->potential);
664 printf("Source\t");
665 printMaxGridMG(ps->source);
666 }
667 /*********************************************************/
printPoissonMG(PoissonMG * ps,int i,int j,int k)668 void printPoissonMG(PoissonMG* ps, int i, int j, int k)
669 {
670 printf("SLICE : %d\t%d\t%d\n",i,j,k);
671 printf("Potential\t");
672 printGridMG(ps->potential,i,j,k);
673 printf("Source\t\t");
674 printGridMG(ps->source,i,j,k);
675 }
676 /*********************************************************/
solveCGPoissonMG(PoissonMG * ps,int max,gdouble acc)677 void solveCGPoissonMG(PoissonMG* ps, int max, gdouble acc)
678 {
679 OperationTypeMG opPotential = getOperationGridMG(ps->potential);
680 OperationTypeMG opSource = getOperationGridMG(ps->source);
681 GridMG* d=getNewGridMG();
682 GridMG* r=getNewGridMG();
683 GridMG* q=getNewGridMG();
684 GridMG* t=getNewGridMG();
685 gdouble deltaNew;
686 gdouble deltaOld;
687 gdouble delta0;
688 gdouble alpha;
689 gdouble beta;
690 gdouble rms;
691 gint i;
692 gdouble scale = 0;
693 gchar tmp[100];
694
695 setOperationGridMG(d,GABEDIT_INTERIOR);
696 setOperationGridMG(r,GABEDIT_INTERIOR);
697 setOperationGridMG(q,GABEDIT_INTERIOR);
698 setOperationGridMG(t,GABEDIT_INTERIOR);
699
700 setOperationGridMG(ps->potential,GABEDIT_INTERIOR);
701 setOperationGridMG(ps->source,GABEDIT_INTERIOR);
702
703 r = residualPoissonMG(ps);
704
705 copyGridMG(d,r);
706 deltaNew = dotGridMG(r,r);
707 delta0 = deltaNew;
708
709 rms = residualNormPoissonMG(ps);
710 tradesBoundaryPoissonMG(ps);
711 copyGridMG(q,r); /* for allocation */
712 copyGridMG(t,r); /* for allocation */
713 scale = (gdouble)1.01/max;
714 progress_orb(0,GABEDIT_PROGORB_COMPMEPGRID,TRUE);
715 for(i=0; i<max && rms>acc ;i++)
716 {
717 if(ps->condition == GABEDIT_CONDITION_PERIODIC) tradesBoundaryPoissonMG(ps);
718 laplacianGridMG(q,d);
719 alpha = deltaNew/dotGridMG(d,q);
720
721 copyGridMG(t,d);
722 multEqualRealGridMG(t,alpha);
723 plusEqualGridMG(ps->potential,t);
724
725 if((i+1)%50 == 0)
726 {
727 r = residualPoissonMG(ps);
728 }
729 else
730 {
731 copyGridMG(t,q);
732 multEqualRealGridMG(t,-alpha);
733 plusEqualGridMG(r,t);
734 }
735 deltaOld = deltaNew;
736 deltaNew = dotGridMG(r,r);
737 beta = deltaNew/deltaOld;
738
739 multEqualRealGridMG(d,beta);
740 plusEqualGridMG(d,r);
741
742 rms = residualNormPoissonMG(ps);
743 /* printf("Solve Poisson by CG i = %d RMS = %f\n",i,rms);*/
744 progress_orb(scale,GABEDIT_PROGORB_COMPMEPGRID,FALSE);
745 sprintf(tmp,_("MEP : Poisson by CG, rms = %f"),rms);
746 setTextInProgress(tmp);
747 if(CancelCalcul)
748 {
749 progress_orb(0,GABEDIT_PROGORB_COMPMEPGRID,TRUE);
750 break;
751 }
752 }
753 setOperationGridMG(ps->potential, opPotential);
754 setOperationGridMG(ps->source, opSource);
755 destroyGridMG(d);
756 destroyGridMG(r);
757 destroyGridMG(q);
758 destroyGridMG(t);
759 }
760 /*********************************************************/
solveMGPoissonMG(PoissonMG * ps,int levelMax)761 gdouble solveMGPoissonMG(PoissonMG* ps, int levelMax)
762 {
763 int level;
764 static int MaxSmmother = 5;
765 gdouble rms = -1;
766 int i;
767
768 if(levelMax<1)
769 {
770 printf(" Error levelMax=%d < 1\n",levelMax);
771 return rms;
772 }
773 smootherPoissonMG(ps,MaxSmmother);
774 if(levelMax==1)
775 {
776 rms = residualNormPoissonMG(ps);
777 return rms;
778 }
779 PoissonMG* poisson=NULL;
780
781 GridMG** e = g_malloc((levelMax-1)*sizeof(GridMG*));
782 GridMG** r = g_malloc((levelMax-1)*sizeof(GridMG*));
783
784
785 level = levelMax -1;
786 i = level - 1;
787 r[i] = residualPoissonMG(ps);
788 setOperationGridMG(r[i], GABEDIT_ALL);
789 restrictionInjectionGridMG(r[i]);
790 e[i] = getNewGridMGFromOldGrid(r[i]);
791 setOperationGridMG(e[i], GABEDIT_ALL);
792 initGridMG(e[i],0.0);
793 poisson = getPoissonMG2(e[i], r[i], ps->condition, NULL);
794 smootherPoissonMG(poisson,MaxSmmother);
795 setOperationGridMG(e[i], GABEDIT_ALL);
796
797
798 for(level = levelMax-2; level>=1; level--)
799 {
800 i = level - 1;
801 r[i] = residualPoissonMG(poisson);
802 setOperationGridMG(r[i], GABEDIT_ALL);
803 restrictionInjectionGridMG(r[i]);
804 e[i] = getNewGridMGFromOldGrid(r[i]);
805 setOperationGridMG(e[i], GABEDIT_ALL);
806 initGridMG(e[i],0.0);
807 poisson = getPoissonMG2(e[i], r[i], ps->condition, NULL);
808 smootherPoissonMG(poisson,MaxSmmother);
809 setOperationGridMG(e[i], GABEDIT_ALL);
810 }
811
812 smootherPoissonMG(poisson,MaxSmmother);
813 setOperationGridMG(e[0], GABEDIT_ALL);
814
815 for(level = 2; level <= levelMax-1; level++)
816 {
817 i = level - 1;
818 prolongationGridMG(e[i-1]);
819 setOperationGridMG(e[i-1],GABEDIT_INTERIOR);
820 setOperationGridMG(e[i],GABEDIT_INTERIOR);
821 plusEqualGridMG(e[i], e[i-1]);
822
823 setOperationGridMG(e[i],GABEDIT_ALL);
824 setOperationGridMG(r[i],GABEDIT_ALL);
825 poisson = getPoissonMG2(e[i], r[i], ps->condition, NULL);
826 smootherPoissonMG(poisson,MaxSmmother);
827 }
828
829
830 i = levelMax-2;
831 /*
832 printf(" avant prolong e[%d] domain\n",i);
833 printDomain(&(e[i]->domain));
834 */
835 prolongationGridMG(e[i]);
836 setOperationGridMG(e[i],GABEDIT_INTERIOR);
837 setOperationGridMG(ps->potential,GABEDIT_INTERIOR);
838 plusEqualGridMG(ps->potential, e[i]);
839 setOperationGridMG(ps->potential,GABEDIT_ALL);
840 smootherPoissonMG(ps,MaxSmmother);
841
842 rms = residualNormPoissonMG(ps);
843 for(i=0;i<levelMax-1;i++)
844 {
845 destroyGridMG(e[i]);
846 destroyGridMG(r[i]);
847 }
848 g_free(r);
849 g_free(e);
850 return rms;
851 }
852 /*********************************************************/
solveMGPoissonMG2(PoissonMG * ps,int levelMax,gdouble acc,int verbose)853 void solveMGPoissonMG2(PoissonMG* ps, int levelMax, gdouble acc, int verbose)
854 {
855 gdouble rms = -1;
856 gint i;
857 if(verbose>=4)
858 {
859 printf("-------------------------------------------\n");
860 printf("Solve Poisson equation by MultiGrid method\n");
861 printf("-------------------------------------------\n");
862 }
863
864 for(i=0;i<500;i++)
865 {
866 rms = solveMGPoissonMG(ps,levelMax);
867 if(verbose>=4)
868 printf("rms = %f\n",rms);
869
870 if(rms>0 && rms<acc)
871 {
872 if(verbose>=4)
873 printf("Number of iterations = %d\n",i+1);
874 break;
875 }
876 }
877 if(verbose>=4)
878 printf("===========================================\n");
879 }
880 /*********************************************************/
solveMGPoissonMG3(PoissonMG * ps,int levelMax,int nIter,gdouble acc,int verbose)881 void solveMGPoissonMG3(PoissonMG* ps, int levelMax, int nIter, gdouble acc, int verbose)
882 {
883 gdouble rms = -1;
884 gint i;
885 gdouble scale = 0;
886 gchar tmp[100];
887 if(verbose>=4)
888 {
889 printf("-------------------------------------------\n");
890 printf("Solve Poisson equation by MultiGrid method\n");
891 printf("-------------------------------------------\n");
892 }
893 progress_orb(0,GABEDIT_PROGORB_COMPGRID,TRUE);
894 setTextInProgress(_("Solve Poisson equation by MultiGrid method, please wait"));
895 scale = (gdouble)1.01/nIter;
896
897
898 for(i=0;i<nIter;i++)
899 {
900 rms = solveMGPoissonMG(ps,levelMax);
901 if(verbose>=4)
902 printf("rms = %f\n",rms);
903
904 if(rms>0 && rms<acc)
905 {
906 if(verbose>=4) printf("Number of iterations = %d\n",i+1);
907 sprintf(tmp,_("MEP : Convergence after %d iterations"),i+1);
908 setTextInProgress(tmp);
909 break;
910 }
911 progress_orb(scale,GABEDIT_PROGORB_COMPMEPGRID,FALSE);
912 sprintf(tmp,_("MEP : Poisson by MultiGrid, rms = %f"),rms);
913 setTextInProgress(tmp);
914 if(CancelCalcul)
915 {
916 progress_orb(0,GABEDIT_PROGORB_COMPMEPGRID,TRUE);
917 break;
918 }
919 }
920 if(verbose>=4)
921 printf("===========================================\n");
922 }
923 /*********************************************************/
solveSmootherPoissonMG2(PoissonMG * ps,int max,int nf)924 void solveSmootherPoissonMG2(PoissonMG* ps, int max, int nf)
925 {
926 printf("In solveSmoother---------------\n");
927 gint n = max/nf;
928 gint i;
929
930 setOperationGridMG(ps->potential,GABEDIT_ALL);
931 setOperationGridMG(ps->source,GABEDIT_ALL);
932 tradesBoundaryPoissonMG(ps);
933
934 for(i=1;i<=n;i++)
935 {
936 smootherPoissonMG(ps,nf);
937 printf(" RMS = %f\n",residualNormPoissonMG(ps));
938 }
939 if(max%nf != 0)
940 smootherPoissonMG(ps,max%nf);
941 }
942 /*********************************************************/
solveSmootherPoissonMG(PoissonMG * ps,int imax,gdouble eps)943 void solveSmootherPoissonMG(PoissonMG* ps, int imax, gdouble eps)
944 {
945 gint i;
946 setOperationGridMG(ps->potential,GABEDIT_ALL);
947 setOperationGridMG(ps->source,GABEDIT_ALL);
948 tradesBoundaryPoissonMG(ps);
949
950 for(i=1;i<=imax;i++)
951 {
952 smootherPoissonMG(ps,1);
953 if( residualNormPoissonMG(ps) <=eps) break;
954 }
955 }
956 /*********************************************************/
957