LORENE
vector_divfree_A_tau.C
1/*
2 * Copyright (c) 2005 Philippe Grandclement
3
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char vector_divfree_A_tau[] = "$Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $" ;
24
25/*
26 * $Id: vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $
27 * $Log: vector_divfree_A_tau.C,v $
28 * Revision 1.4 2014/10/13 08:53:45 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.3 2014/10/06 15:13:20 j_novak
32 * Modified #include directives to use c++ syntax.
33 *
34 * Revision 1.2 2013/06/05 15:10:43 j_novak
35 * Suppression of FINJAC sampling in r. This Jacobi(0,2) base is now
36 * available by setting colloc_r to BASE_JAC02 in the Mg3d constructor.
37 *
38 * Revision 1.1 2008/08/27 09:01:27 jl_cornou
39 * Methods for solving Dirac systems for divergence free vectors
40 *
41 * Revision 1.6 2007/12/14 10:19:34 jl_cornou
42 * *** empty log message ***
43 *
44 * Revision 1.4 2005/11/24 14:07:54 j_novak
45 * Use of Matrice::annule_hard()
46 *
47 * Revision 1.3 2005/08/26 14:02:41 p_grandclement
48 * Modification of the elliptic solver that matches with an oscillatory exterior solution
49 * small correction in Poisson tau also...
50 *
51 * Revision 1.2 2005/08/25 12:16:01 p_grandclement
52 * *** empty log message ***
53 *
54 *
55 * $Header: /cvsroot/Lorene/C++/Source/Tensor/vector_divfree_A_tau.C,v 1.4 2014/10/13 08:53:45 j_novak Exp $
56 *
57 */
58
59// C headers
60#include <cassert>
61#include <cmath>
62
63// Lorene headers
64#include "tensor.h"
65#include "diff.h"
66#include "proto.h"
67#include "param.h"
68#include "matrice.h"
69#include "type_parite.h"
70
71
72namespace Lorene {
74
75 const Map_af* mp_aff = dynamic_cast<const Map_af*>(mp) ;
76 assert(mp_aff != 0x0) ; // Only affine mapping for the moment
77
78 const Mg3d& mgrid = *mp_aff->get_mg();
79 assert((mgrid.get_type_r(0) == RARE) || (mgrid.get_type_r(0) == FIN));
80 if (aaa.get_etat() == ETATZERO) {
81 eta_tilde = 0;
82 vr = 0 ;
83 return ;
84 }
85 assert(aaa.get_etat() != ETATNONDEF) ;
86 int nz = mgrid.get_nzone();
87 int nzm1 = nz - 1 ;
88 bool ced = (mgrid.get_type_r(nzm1) == UNSURR) ;
89 int n_shell = ced ? nz-2 : nzm1 ;
90 int nz_bc = nzm1 ;
91 if (par_bc != 0x0)
92 if (par_bc->get_n_int() > 0) nz_bc = par_bc->get_int();
94 //bool cedbc (ced && (nz_bc == nzm1));
95
96//#ifndef NDEBUG
97// if (!cedbc) {
98// assert(par_bc != 0x0) ;
99// assert(par_bc->get_n_tbl_mod() >= 3) ;
100// }
101//#endif
102 int nt = mgrid.get_nt(0) ;
103 int np = mgrid.get_np(0) ;
104 int nr ;
105
106 Scalar source = aaa ;
108 source_coq.annule_domain(0);
109
110 int dzpuis = source.get_dzpuis();
111
112 if (ced) source_coq.annule_domain(nzm1) ;
113 source_coq.mult_r() ;
114 source.set_spectral_va().ylm() ;
115 source_coq.set_spectral_va().ylm() ;
116 Base_val base = source.get_spectral_base() ;
117 base.mult_x() ;
118
119 eta_tilde.annule_hard() ;
120 eta_tilde.set_spectral_base(base) ;
121 eta_tilde.set_spectral_va().set_etat_cf_qcq() ;
122 eta_tilde.set_spectral_va().c_cf->annule_hard() ;
123 vr.annule_hard() ;
124 vr.set_spectral_base(base) ;
125 vr.set_spectral_va().set_etat_cf_qcq() ;
126 vr.set_spectral_va().c_cf->annule_hard() ;
127 Mtbl_cf& meta = *eta_tilde.set_spectral_va().c_cf ;
128 Mtbl_cf& mvr = *vr.set_spectral_va().c_cf ;
129
130
131 int base_r ;
132
133
134 // Determination of the size of the systeme :
135 int size = 0 ;
136 int max_nr = 0 ;
137 for (int l=0 ; l<nz ; l++) {
138 nr = mgrid.get_nr(l) ;
139 size += 2*nr ;
140 if (nr > max_nr)
141 max_nr = nr ;
142 }
143
144 Matrice systeme (size, size) ;
145 systeme.set_etat_qcq() ;
146 Tbl sec_membre (size) ;
147
148 np = mgrid.get_np(0) ;
149 nt = mgrid.get_nt(0) ;
150
151
152 double alpha, beta ;
153 int l_quant, m_quant ;
154
155 // On bosse pour chaque l, m :
156 for (int k=0 ; k<np+1 ; k++)
157 for (int j=0 ; j<nt ; j++)
158 if (nullite_plm(j, nt, k, np, base) == 1) {
159
160 systeme.annule_hard() ;
161 sec_membre.annule_hard() ;
162
163 int column_courant = 0 ;
164 int ligne_courant = 0 ;
165
166 //--------------------------
167 // NUCLEUS
168 //--------------------------
169
170 nr = mgrid.get_nr(0) ;
171 alpha = (*mp_aff).get_alpha()[0] ;
172 base.give_quant_numbers (0, k, j, m_quant, l_quant, base_r) ;
173 Diff_dsdx odn(base_r, nr) ; const Matrice& mdn = odn.get_matrice() ;
174 Diff_sx osn(base_r, nr) ; const Matrice& msn = osn.get_matrice() ;
175
176
177
178 int nbr_cl = 0 ;
179 // RARE case
180 if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
181 // regularity conditions for eta :
182 if (l_quant > 1) {
183 nbr_cl = 1 ;
184 if (l_quant%2==0) {
185 //Even case
186 for (int col=0 ; col<nr ; col++) {
187 if (col%2==0) {
189 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
190 else {
192 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
193 }
194 }
195 else {
196 //Odd case
197 for (int col=0 ; col<nr ; col++) {
198 if (col%2==0) {
200 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
201 else {
203 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
204 }
205 }
206 }
207 }
208
209 // R_JACO02 case
210 else {
211 assert (base_r == R_JACO02) ;
212 // regularity conditions for eta :
213 if (l_quant == 0) {
214 nbr_cl = 1 ;
215 for (int col=0 ; col<nr ; col++) {
216 systeme.set(ligne_courant, col+column_courant) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1);
218 }
219 }
220 else if (l_quant == 1) {
221 nbr_cl = 1 ;
222 for (int col=0 ; col<nr ; col++) {
223 systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
225 }
226 }
227 else {
228 nbr_cl = 2 ;
229 for (int col=0 ; col<nr ; col++) {
230 systeme.set(ligne_courant, col+column_courant) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
232 systeme.set(ligne_courant+1, col+column_courant) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
233 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
234 }
235 }
236 }
238
239 // Premiere partie de l'operateur :
240 for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
241 for (int col=0 ; col<nr ; col++) {
243 sec_membre.set(lig+ligne_courant) = alpha*alpha*(*source.get_spectral_va().c_cf)(0, k, j, lig) ;
245 // systeme.set(lig+ligne_courant+nr) = 0 ;
246 }
247 }
248
249
250 ligne_courant += nr-1-nbr_cl ;
251
252
253 // RARE case
254 if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
255 // regularity conditions for vr :
256 if (l_quant > 1) {
257 nbr_cl = 1 ;
258 if (l_quant%2==0) {
259 //Even case
260 for (int col=0 ; col<nr ; col++) {
261 if (col%2==0) {
263 systeme.set(ligne_courant, col+column_courant+nr) = 1 ; }
264 else {
266 systeme.set(ligne_courant, col+column_courant+nr) = -1 ; }
267 }
268 }
269
270 else {
271 //Odd case
272 for (int col=0 ; col<nr ; col++) {
273 if (col%2==0) {
275 systeme.set(ligne_courant, col+column_courant+nr) = 2*col+1 ; }
276 else {
278 systeme.set(ligne_courant, col+column_courant+nr) = -(2*col+1) ; }
279 }
280 }
281 }
282 }
283
284 // R_JACO02 case
285 else {
286 assert (base_r == R_JACO02) ;
287 // regularity conditions for vr :
288 if (l_quant == 0) {
289 nbr_cl = 1 ;
290 for (int col=0 ; col<nr ; col++) {
292 systeme.set(ligne_courant, col+column_courant+nr) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
293 }
294 }
295 else if (l_quant == 1) {
296 nbr_cl = 1 ;
297 for (int col=0 ; col<nr ; col++) {
299 systeme.set(ligne_courant, col+column_courant+nr) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
300 }
301 }
302 else {
303 nbr_cl = 2 ;
304 for (int col=0 ; col<nr ; col++) {
306 systeme.set(ligne_courant, col+column_courant+nr) = (col+1)*(col+2)/double(2)*(1-2*(col%2)) ;
308 systeme.set(ligne_courant+1, col+column_courant+nr) = col*(col+1)*(col+2)*(col+3)/double(12)*(2*(col%2)-1) ;
309 }
310 }
311 }
313
314 // Deuxieme partie de l'operateur
315 for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
316 for (int col=0 ; col<nr ; col++) {
317 systeme.set(lig+ligne_courant,col+column_courant) = - l_quant*(l_quant+1)*msn(lig,col) ;
318 sec_membre.set(lig+ligne_courant) = 0 ;
320 // systeme.set(lig+ligne_courant+nr) = 0 ;
321 }
322 }
323
324
325 ligne_courant += nr-1-nbr_cl ;
326
327
328 // Le raccord pour eta
329 for (int col=0 ; col<nr ; col++) {
330 if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
331 // La fonction eta
334 // Sa derivee :
335 if (l_quant%2==0) {
336 systeme.set(ligne_courant+1, col+column_courant) = 4*col*col/alpha ;
337 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
338 }
339 else {
340 systeme.set(ligne_courant+1, col+column_courant) = (2*col+1)*(2*col+1)/alpha ;
341 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
342 }
343 }
344 else { // Cas R_JACO02
345 assert( base_r == R_JACO02 ) ;
346 // La fonction eta
349 // Sa derivee :
350 systeme.set(ligne_courant+1, col+column_courant) = col*(col+3)/double(2)/alpha ;
351 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
352 }
353 }
354
355 ligne_courant += 2 ;
356 // Le raccord pour vr
357 for (int col=0 ; col<nr ; col++) {
358 if (source.get_mp().get_mg()->get_type_r(0) == RARE) {
359 // La fonction vr
362 // Sa derivee :
363 if (l_quant%2==0) {
365 systeme.set(ligne_courant+1, col+column_courant+nr) = 4*col*col/alpha ;
366 }
367 else {
369 systeme.set(ligne_courant+1, col+column_courant+nr) = (2*col+1)*(2*col+1)/alpha ;
370 }
371 }
372 else { // Cas R_JACO02
373 assert( base_r == R_JACO02 ) ;
374 // La fonction vr
377 // Sa derivee :
379 systeme.set(ligne_courant+1, col+column_courant+nr) = col*(col+3)/double(2)/alpha ;
380 }
381 }
382
383 ligne_courant -= 2 ; // On retourne pour le raccord dans la prochaine zone
384
385 column_courant += 2*nr ; // On va changer de zone
386
387 //--------------------------
388 // SHELLS
389 //--------------------------
390 for (int l=1 ; l<nz-1 ; l++) {
391
392 nr = mgrid.get_nr(l) ;
393 alpha = (*mp_aff).get_alpha()[l] ;
394 beta = (*mp_aff).get_beta()[l] ;
395
396 base.give_quant_numbers (l, k, j, m_quant, l_quant, base_r) ;
397 Diff_id ois(base_r, nr) ; const Matrice& mis = ois.get_matrice() ;
398 Diff_xdsdx oxds(base_r, nr) ; const Matrice& mxds = oxds.get_matrice() ;
399
400 // matching with previous domain :
401 for (int col=0 ; col<nr ; col++) {
402 // La fonction eta
403 if (col%2==0) {
405 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
406 else {
408 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
409 // Sa derivee :
410 if (col%2==0) {
412 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
413 else {
414 systeme.set(ligne_courant+1, col+column_courant) = -col*col/alpha ;
415 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
416 }
417 ligne_courant += 2 ;
418
419 // matching with previous domain :
420 for (int col=0 ; col<nr ; col++) {
421 // La fonction vr
422 if (col%2==0) {
424 systeme.set(ligne_courant, col+column_courant+nr) = -1 ; }
425 else {
427 systeme.set(ligne_courant, col+column_courant+nr) = 1 ; }
428 // Sa derivee :
429 if (col%2==0) {
431 systeme.set(ligne_courant+1, col+column_courant+nr) = col*col/alpha ; }
432 else {
434 systeme.set(ligne_courant+1, col+column_courant+nr) = -col*col/alpha ; }
435 }
436 ligne_courant += 2 ;
437
438
439 // L'operateur (premiere et deuxieme parties) :
440
441 // source must be multiplied by (x+echelle)^2
442 Tbl source_aux(nr) ;
443 source_aux.set_etat_qcq() ;
444 for (int i=0 ; i<nr ; i++)
445 source_aux.set(i) = (*source.get_spectral_va().c_cf)(l,k,j,i)*alpha*alpha ;
448 multx_1d(nr, &xso.t, R_CHEB) ;
449 multx_1d(nr, &xxso.t, R_CHEB) ;
450 multx_1d(nr, &xxso.t, R_CHEB) ;
451 source_aux = beta*beta/alpha/alpha*source_aux+2*beta/alpha*xso+xxso ;
452
453 for (int lig=0 ; lig<nr-2 ; lig++) {
454 for (int col=0 ; col<nr ; col++) {
458 systeme.set(lig+ligne_courant+nr-2, col+column_courant) = -l_quant*(l_quant+1) ;
460 sec_membre.set(lig+ligne_courant+nr-2) = 0 ; }
461 }
462
463
464 ligne_courant += 2*nr-4 ;
465 // Matching with the next domain :
466 for (int col=0 ; col<nr ; col++) {
467 // La fonction eta
470 // Sa derivee :
472 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ;
473 // La fonction vr
475 systeme.set(ligne_courant+2, col+column_courant+nr) = 1 ;
476 // Sa derivee
478 systeme.set(ligne_courant+3, col+column_courant+nr) = col*col/alpha ;
479 }
480
481 column_courant += 2*nr ;
482 }
483
484
485 //--------------------------
486 // ZEC
487 //--------------------------
488 nr = mgrid.get_nr(nz-1) ;
489 alpha = (*mp_aff).get_alpha()[nz-1] ;
490 beta = (*mp_aff).get_beta()[nz-1] ;
491
492 base.give_quant_numbers (nz-1, k, j, m_quant, l_quant, base_r) ;
493 //work = new Matrice(laplacien_mat(nr, l_quant, 0., dzpuis, base_r)) ;
494 Diff_dsdx odzec(base_r, nr) ; const Matrice& mdzec = odzec.get_matrice() ;
495 Diff_sx oszec(base_r, nr) ; const Matrice& mszec = oszec.get_matrice() ;
496
497
498 // Matching with the previous domain :
499 for (int col=0 ; col<nr ; col++) {
500 // La fonction eta
501 if (col%2==0) {
503 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
504 else {
506 systeme.set(ligne_courant, col+column_courant+nr) = 0 ; }
507 // Sa derivee :
508 if (col%2==0) {
509 systeme.set(ligne_courant+1, col+column_courant) = -4*alpha*col*col ;
510 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
511 else {
512 systeme.set(ligne_courant+1, col+column_courant) = 4*alpha*col*col ;
513 systeme.set(ligne_courant+1, col+column_courant+nr) = 0 ; }
514 // La fonction vr
515 if (col%2==0) {
517 systeme.set(ligne_courant+2, col+column_courant+nr) = 0 ; }
518 else {
520 systeme.set(ligne_courant+2, col+column_courant+nr) = 0 ; }
521 // Sa derivee :
522 if (col%2==0) {
523 systeme.set(ligne_courant+3, col+column_courant) = -4*alpha*col*col ;
524 systeme.set(ligne_courant+3, col+column_courant+nr) = 0 ; }
525 else {
526 systeme.set(ligne_courant+3, col+column_courant) = 4*alpha*col*col ;
527 systeme.set(ligne_courant+3, col+column_courant+nr) = 0 ; }
528 }
529 ligne_courant += 4 ;
530
531 // Regularity and BC at infinity ?
532 nbr_cl =0 ;
533 switch (dzpuis) {
534 case 4 :
535 if (l_quant==0) {
536 nbr_cl = 1 ;
537 // Only BC at infinity :
538 for (int col=0 ; col<nr ; col++) {
542 systeme.set(ligne_courant+1, col+column_courant+nr) = 1 ; }
543 }
544 else {
545 nbr_cl = 2 ;
546 // BC at infinity :
547 for (int col=0 ; col<nr ; col++) {
551 systeme.set(ligne_courant+1, col+column_courant+nr) = 1; }
552 // Regularity :
553 for (int col=0 ; col<nr ; col++) {
554 systeme.set(ligne_courant+2, col+column_courant) = -4*alpha*col*col ;
555 systeme.set(ligne_courant+2, col+column_courant+nr) = 0 ;
557 systeme.set(ligne_courant+3, col+column_courant+nr) = -4*alpha*col*col ; }
558 }
559 break ;
560
561 case 3 :
562 nbr_cl = 1 ;
563 // Only BC at infinity :
564 for (int col=0 ; col<nr ; col++) {
568 systeme.set(ligne_courant+1, col+column_courant+nr) = 1 ; }
569 break ;
570
571 case 2 :
572 if (l_quant==0) {
573 nbr_cl = 1 ;
574 // Only BC at infinity :
575 for (int col=0 ; col<nr ; col++) {
579 systeme.set(ligne_courant+1, col+column_courant+1) = 1 ; }
580 }
581 break ;
582 default :
583 cout << "Unknown dzpuis in vector_divfree_A ..." << endl ;
584 abort() ;
585 }
586
588
589 // Multiplication of the source :
590 double indic = 1 ;
591 switch (dzpuis) {
592 case 4 :
593 indic = alpha*alpha ;
594 break ;
595 case 3 :
596 indic = alpha ;
597 break ;
598 default :
599 break ;
600 }
601
602 // L'operateur :
603 for (int lig=0 ; lig<nr-1-nbr_cl ; lig++) {
604 for (int col=0 ; col<nr ; col++) {
607 sec_membre.set(lig+ligne_courant) = indic*(*source.get_spectral_va().c_cf)(nz-1, k, j, lig) ;
608 systeme.set(lig+ligne_courant+nr-1-nbr_cl,col+column_courant)= - l_quant*(l_quant+1)*mszec(lig,col) ;
610 sec_membre.set(lig+ligne_courant+nr-1-nbr_cl) = 0 ; }
611 }
612
613
614 // Solving the system:
615 systeme.set_band (max_nr, max_nr) ;
616 systeme.set_lu() ;
617 Tbl soluce (systeme.inverse(sec_membre)) ;
618
619 // On range :
620 int conte = 0 ;
621 for (int l=0 ; l<nz ; l++) {
622 nr = mgrid.get_nr(l) ;
623 for (int i=0 ; i<nr ; i++) {
624 meta.set(l,k,j,i) = soluce(conte) ;
625 mvr.set(l,k,i,j) = soluce(conte+nr) ;
626 conte ++ ;
627 }
628 }
629}
630if (eta_tilde.set_spectral_va().c != 0x0)
631 delete eta_tilde.set_spectral_va().c ;
632 eta_tilde.set_spectral_va().c = 0x0 ;
633 eta_tilde.set_spectral_va().ylm_i() ;
634
635 if (vr.set_spectral_va().c != 0x0)
636 delete vr.set_spectral_va().c ;
637 vr.set_spectral_va().c = 0x0 ;
638 vr.set_spectral_va().ylm_i() ;
639}
640}
Bases of the spectral expansions.
Definition base_val.h:322
void mult_x()
The basis is transformed as with a multiplication by .
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:129
Class for the elementary differential operator Identity (see the base class Diff ).
Definition diff.h:210
Class for the elementary differential operator division by (see the base class Diff ).
Definition diff.h:329
Class for the elementary differential operator (see the base class Diff ).
Definition diff.h:409
Time evolution with partial storage (*** under development ***).
Definition evolution.h:371
Affine radial mapping.
Definition map.h:2027
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Matrix handling.
Definition matrice.h:152
Multi-domain grid.
Definition grilles.h:273
Coefficients storage for the multi-domain spectral method.
Definition mtbl_cf.h:186
Parameter storage.
Definition param.h:125
Tensor field of valence 0 (or component of a tensorial field).
Definition scalar.h:387
Basic array class.
Definition tbl.h:161
void sol_Dirac_A_tau(const Scalar &aaa, Scalar &eta, Scalar &vr, const Param *par_bc=0x0) const
Solves via a tau method a system of two-coupled first-order PDEs obtained from the divergence-free co...
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
#define R_CHEB
base de Chebychev ordinaire (fin)
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
Definition tensor.h:295
Lorene prototypes.
Definition app_hor.h:64