LORENE
eos_bifluid.C
1/*
2 * Methods of the class Eos_bifluid.
3 *
4 * (see file eos_bifluid.h for documentation).
5 */
6
7/*
8 * Copyright (c) 2001 Jerome Novak
9 *
10 * This file is part of LORENE.
11 *
12 * LORENE is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28
29char eos_bifluid_C[] = "$Header: /cvsroot/Lorene/C++/Source/Eos/eos_bifluid.C,v 1.22 2015/06/12 12:38:24 j_novak Exp $" ;
30
31/*
32 * $Id: eos_bifluid.C,v 1.22 2015/06/12 12:38:24 j_novak Exp $
33 * $Log: eos_bifluid.C,v $
34 * Revision 1.22 2015/06/12 12:38:24 j_novak
35 * Implementation of the corrected formula for the quadrupole momentum.
36 *
37 * Revision 1.21 2015/06/11 14:41:59 a_sourie
38 * Corrected minor bug
39 *
40 * Revision 1.20 2015/06/11 13:50:19 j_novak
41 * Minor corrections
42 *
43 * Revision 1.19 2015/06/10 14:39:17 a_sourie
44 * New class Eos_bf_tabul for tabulated 2-fluid EoSs and associated functions for the computation of rotating stars with such EoSs.
45 *
46 * Revision 1.18 2014/10/13 08:52:52 j_novak
47 * Lorene classes and functions now belong to the namespace Lorene.
48 *
49 * Revision 1.17 2014/04/25 10:43:51 j_novak
50 * The member 'name' is of type string now. Correction of a few const-related issues.
51 *
52 * Revision 1.16 2008/08/19 06:42:00 j_novak
53 * Minor modifications to avoid warnings with gcc 4.3. Most of them concern
54 * cast-type operations, and constant strings that must be defined as const char*
55 *
56 * Revision 1.15 2006/03/10 08:38:55 j_novak
57 * Use of C++-style casts.
58 *
59 * Revision 1.14 2004/09/01 16:12:30 r_prix
60 * (hopefully) fixed seg-fault bug with my inconsistent treatment of eos-bifluid 'name'
61 * (was char-array, now char*)
62 *
63 * Revision 1.13 2004/09/01 09:50:34 r_prix
64 * adapted to change in read_variable() for strings
65 *
66 * Revision 1.12 2003/12/17 23:12:32 r_prix
67 * replaced use of C++ <string> by standard ANSI char* to be backwards compatible
68 * with broken compilers like MIPSpro Compiler 7.2 on SGI Origin200. ;-)
69 *
70 * Revision 1.11 2003/12/05 15:09:47 r_prix
71 * adapted Eos_bifluid class and subclasses to use read_variable() for
72 * (formatted) file-reading.
73 *
74 * Revision 1.10 2003/12/04 14:17:26 r_prix
75 * new 2-fluid EOS subtype 'typeos=5': this is identical to typeos=0
76 * (analytic EOS), but we perform the EOS inversion "slow-rot-style",
77 * i.e. we don't switch to a 1-fluid EOS when one fluid vanishes.
78 *
79 * Revision 1.9 2003/11/18 18:28:38 r_prix
80 * moved particle-masses m_1, m_2 of the two fluids into class eos_bifluid (from eos_bf_poly)
81 *
82 * Revision 1.8 2003/10/03 15:58:46 j_novak
83 * Cleaning of some headers
84 *
85 * Revision 1.7 2002/10/16 14:36:35 j_novak
86 * Reorganization of #include instructions of standard C++, in order to
87 * use experimental version 3 of gcc.
88 *
89 * Revision 1.6 2002/05/31 16:13:36 j_novak
90 * better inversion for eos_bifluid
91 *
92 * Revision 1.5 2002/05/10 09:55:27 j_novak
93 * *** empty log message ***
94 *
95 * Revision 1.4 2002/05/10 09:26:52 j_novak
96 * Added new class Et_rot_mag for magnetized rotating neutron stars (under development)
97 *
98 * Revision 1.3 2002/01/03 15:30:27 j_novak
99 * Some comments modified.
100 *
101 * Revision 1.2 2001/12/04 21:27:53 e_gourgoulhon
102 *
103 * All writing/reading to a binary file are now performed according to
104 * the big endian convention, whatever the system is big endian or
105 * small endian, thanks to the functions fwrite_be and fread_be
106 *
107 * Revision 1.1.1.1 2001/11/20 15:19:27 e_gourgoulhon
108 * LORENE
109 *
110 * Revision 1.5 2001/10/10 13:49:53 eric
111 * Modif Joachim: &(Eos_bifluid::...) --> &Eos_bifluid::...
112 * pour conformite au compilateur HP.
113 *
114 * Revision 1.4 2001/08/31 15:48:11 novak
115 * The flag tronc has been added to ar_ent... functions
116 *
117 * Revision 1.3 2001/08/27 12:23:40 novak
118 * The Cmp arguments delta2 put to const
119 *
120 * Revision 1.2 2001/08/27 09:52:49 novak
121 * Use of new variable delta2
122 *
123 * Revision 1.1 2001/06/21 15:21:47 novak
124 * Initial revision
125 *
126 *
127 * $Header: /cvsroot/Lorene/C++/Source/Eos/eos_bifluid.C,v 1.22 2015/06/12 12:38:24 j_novak Exp $
128 *
129 */
130
131// Headers C
132#include <cstdlib>
133#include <cmath>
134
135// Headers Lorene
136#include "eos_bifluid.h"
137#include "cmp.h"
138#include "utilitaires.h"
139
140 //--------------//
141 // Constructors //
142 //--------------//
143
144
145// Standard constructor without name
146// ---------------------------------
147namespace Lorene {
149 name("EoS bi-fluid"), m_1(1), m_2(1)
150{ }
151
152// Standard constructor with name and masses
153// ---------------------------------
154Eos_bifluid::Eos_bifluid(const char* name_i, double mass1, double mass2) :
155 name(name_i), m_1(mass1), m_2(mass2)
156{ }
157
158// Copy constructor
159// ----------------
161 name(eos_i.name), m_1(eos_i.m_1), m_2(eos_i.m_2)
162{ }
163
164// Constructor from a binary file
165// ------------------------------
167{
168 char dummy [MAX_EOSNAME];
169 if (fread(dummy, sizeof(char),MAX_EOSNAME, fich) == 0)
170 cerr << "Reading problem in " << __FILE__ << endl ;
171 name = dummy ;
172 fread_be(&m_1, sizeof(double), 1, fich) ;
173 fread_be(&m_2, sizeof(double), 1, fich) ;
174
175}
176
177// Constructor from a formatted file
178// ---------------------------------
179Eos_bifluid::Eos_bifluid(const char *fname)
180{
181 int res=0;
182 char* dummy = 0x0 ;
183 char* char_name = const_cast<char*>("name");
184 char* char_m1 = const_cast<char*>("m_1") ;
185 char* char_m2 = const_cast<char*>("m_2") ;
186 res += read_variable (fname, char_name, &dummy);
187 res += read_variable (fname, char_m1, m_1);
188 res += read_variable (fname, char_m2, m_2);
189
190 name = dummy ;
191
192 free(dummy) ;
193
194 if (res)
195 {
196 cerr << "ERROR: Eos_bifluid(const char*): could not read either of \
197the variables 'name', 'm_1, or 'm_2' from file " << fname << endl;
198 exit (-1);
199 }
200
201}
202
203// Constructor from a formatted file
204// ---------------------------------
206
207 string aname ;
208
209 // EOS identificator :
210 char blabla[80] ;
211 fich >> aname;
212 name = aname ;
213 fich.getline(blabla, 80) ;
214
215 fich >> m_1 ; fich.getline(blabla, 80) ;
216 fich >> m_2 ; fich.getline(blabla, 80) ;
217}
218
219
220
221
222 //--------------//
223 // Destructor //
224 //--------------//
225
228
229 //--------------//
230 // Assignment //
231 //--------------//
233
234 name = eosi.name ;
235 m_1 = eosi.m_1;
236 m_2 = eosi.m_2;
237
238}
239
240
241 //------------//
242 // Outputs //
243 //------------//
244
245void Eos_bifluid::sauve(FILE* fich) const
246{
247 assert(name.size() < MAX_EOSNAME) ;
248 char dummy [MAX_EOSNAME];
249 int ident = identify() ;
250
251 fwrite_be(&ident, sizeof(int), 1, fich) ;
252
253 strncpy (dummy, name.c_str(), MAX_EOSNAME);
254 dummy[MAX_EOSNAME-1] = 0;
255 if (fwrite(dummy, sizeof(char), MAX_EOSNAME, fich ) == 0)
256 cerr << "Writing problem in " << __FILE__ << endl ;
257
258 fwrite_be(&m_1, sizeof(double), 1, fich) ;
259 fwrite_be(&m_2, sizeof(double), 1, fich) ;
260
261}
262
263
264ostream& operator<<(ostream& ost, const Eos_bifluid& eqetat) {
265 ost << eqetat.get_name() << endl ;
266 ost << " Mean particle 1 mass : " << eqetat.get_m1() << " m_B" << endl ;
267 ost << " Mean particle 2 mass : " << eqetat.get_m2() << " m_B" << endl ;
268
269 eqetat >> ost ;
270 return ost ;
271}
272
273
274 //-------------------------------//
275 // Computational routines //
276 //-------------------------------//
277
278// Complete computational routine giving all thermo variables
279//-----------------------------------------------------------
280
281void Eos_bifluid::calcule_tout(const Cmp& ent1, const Cmp& ent2,
282 const Cmp& delta2, Cmp& nbar1, Cmp& nbar2,
283 Cmp& ener, Cmp& press,
284 int nzet, int l_min) const {
285
286 assert(ent1.get_etat() != ETATNONDEF) ;
287 assert(ent2.get_etat() != ETATNONDEF) ;
288 assert(delta2.get_etat() != ETATNONDEF) ;
289
290 const Map* mp = ent1.get_mp() ; // Mapping
291 assert(mp == ent2.get_mp()) ;
292 assert(mp == delta2.get_mp()) ;
293 assert(mp == ener.get_mp()) ;
294
295 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
296 nbar1.set_etat_zero() ;
297 nbar2.set_etat_zero() ;
298 ener.set_etat_zero() ;
299 press.set_etat_zero() ;
300 return ;
301 }
302 nbar1.allocate_all() ;
303 nbar2.allocate_all() ;
304 ener.allocate_all() ;
305 press.allocate_all() ;
306
307 const Mg3d* mg = mp->get_mg() ; // Multi-grid
308
309 int nz = mg->get_nzone() ; // total number of domains
310
311 // resu is set to zero in the other domains :
312
313 if (l_min > 0) {
314 nbar1.annule(0, l_min-1) ;
315 nbar2.annule(0, l_min-1) ;
316 ener.annule(0, l_min-1) ;
317 press.annule(0, l_min-1) ;
318 }
319
320 if (l_min + nzet < nz) {
321 nbar1.annule(l_min + nzet, nz - 1) ;
322 nbar2.annule(l_min + nzet, nz - 1) ;
323 ener.annule(l_min + nzet, nz - 1) ;
324 press.annule(l_min + nzet, nz - 1) ;
325 }
326
327 int np0 = mp->get_mg()->get_np(0) ;
328 int nt0 = mp->get_mg()->get_nt(0) ;
329 for (int l=l_min; l<l_min+nzet; l++) {
330 assert(mp->get_mg()->get_np(l) == np0) ;
331 assert(mp->get_mg()->get_nt(l) == nt0) ;
332 }
333
334 //**********************************************************************
335 //RP: for comparison with slow-rotation, we might have to treat the
336 // 1-fluid region somewhat differently...
337 bool slow_rot_style = false; // off by default
338
339 if ( identify() == 2 ) // only applies if newtonian 2-fluid polytrope
340 {
341 const Eos_bf_poly_newt *this_eos = dynamic_cast<const Eos_bf_poly_newt*>(this);
342 if (this_eos -> get_typeos() == 5)
343 {
344 slow_rot_style = true;
345 cout << "DEBUG: using EOS-inversion slow-rot-style!! " << endl;
346 }
347 }
348
349 //**********************************************************************
350
351 for (int k=0; k<np0; k++) {
352 for (int j=0; j<nt0; j++) {
353 bool inside = true ;
354 bool inside1 = true ;
355 bool inside2 = true ;
356 for (int l=l_min; l< l_min + nzet; l++) {
357 for (int i=0; i<mp->get_mg()->get_nr(l); i++) {
358 double xx, xx1, xx2, n1, n2 ;
359 xx1 = ent1(l,k,j,i) ;
360 xx2 = ent2(l,k,j,i) ;
361 xx = delta2(l,k,j,i) ;
362 if (inside) {
363 inside = (!nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
364
365 // inside1 = ((n1 > 0.)&&(xx1>0.)) ;
366 inside1 = (n1 > 0.) ;
367 // inside2 = ((n2 > 0.)&&(xx2>0.)) ;
368 inside2 = (n2 > 0.) ;
369
370 // slowrot special treatment follows here.
371 if (slow_rot_style)
372 {
373 inside = true; // no 1-fluid transition!
374 n1 = (n1 > 0) ? n1: 0; // make sure only positive densities
375 n2 = (n2 > 0) ? n2: 0;
376 }
377 }
378 if (inside) {
379 nbar1.set(l,k,j,i) = n1 ;
380 nbar2.set(l,k,j,i) = n2 ;
381 }
382 else {
383 if (inside1) {
384 n1 = nbar_ent_p1(xx1) ;
385 inside1 = (n1 > 0.) ;
386 }
387 if (inside2) {
388 n2 = nbar_ent_p2(xx2) ;
389 inside2 = (n2 > 0.) ;
390 }
391 if (!inside1) n1 = 0. ;
392 if (!inside2) n2 = 0. ;
393 nbar1.set(l,k,j,i) = n1 ;
394 nbar2.set(l,k,j,i) = n2 ;
395 }
396
397 ener.set(l,k,j,i) = ener_nbar_p(n1, n2, xx) ;
398 press.set(l,k,j,i) = press_nbar_p(n1, n2, xx) ;
399
400 }
401 }
402 }
403
404 }
405}
406
407
408// Baryon density from enthalpy
409//------------------------------
410
411void Eos_bifluid::nbar_ent(const Cmp& ent1, const Cmp& ent2, const Cmp& delta2,
412 Cmp& nbar1, Cmp& nbar2, int nzet, int l_min) const {
413
414 assert(ent1.get_etat() != ETATNONDEF) ;
415 assert(ent2.get_etat() != ETATNONDEF) ;
416 assert(delta2.get_etat() != ETATNONDEF) ;
417
418 const Map* mp = ent1.get_mp() ; // Mapping
419 assert(mp == ent2.get_mp()) ;
420 assert(mp == delta2.get_mp()) ;
421 assert(mp == nbar1.get_mp()) ;
422
423 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
424 nbar1.set_etat_zero() ;
425 nbar2.set_etat_zero() ;
426 return ;
427 }
428 nbar1.allocate_all() ;
429 nbar2.allocate_all() ;
430
431 const Mg3d* mg = mp->get_mg() ; // Multi-grid
432
433 int nz = mg->get_nzone() ; // total number of domains
434
435 // resu is set to zero in the other domains :
436
437 if (l_min > 0) {
438 nbar1.annule(0, l_min-1) ;
439 nbar2.annule(0, l_min-1) ;
440 }
441
442 if (l_min + nzet < nz) {
443 nbar1.annule(l_min + nzet, nz - 1) ;
444 nbar2.annule(l_min + nzet, nz - 1) ;
445 }
446
447 int np0 = mp->get_mg()->get_np(0) ;
448 int nt0 = mp->get_mg()->get_nt(0) ;
449 for (int l=l_min; l<l_min+nzet; l++) {
450 assert(mp->get_mg()->get_np(l) == np0) ;
451 assert(mp->get_mg()->get_nt(l) == nt0) ;
452 }
453
454 for (int k=0; k<np0; k++) {
455 for (int j=0; j<nt0; j++) {
456 bool inside = true ;
457 bool inside1 = true ;
458 bool inside2 = true ;
459 for (int l=l_min; l< l_min + nzet; l++) {
460 for (int i=0; i<mp->get_mg()->get_nr(l); i++) {
461 double xx, xx1, xx2, n1, n2 ;
462 xx1 = ent1(l,k,j,i) ;
463 xx2 = ent2(l,k,j,i) ;
464 xx = delta2(l,k,j,i) ;
465 if (inside) {
466 inside = (!nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
467 inside1 = ((n1 > 0.)&&(xx1>0.)) ;
468 inside2 = ((n2 > 0.)&&(xx2>0.)) ;
469 }
470 if (inside) {
471 nbar1.set(l,k,j,i) = n1 ;
472 nbar2.set(l,k,j,i) = n2 ;
473 }
474 else {
475 if (inside1) {
476 n1 = nbar_ent_p1(xx1) ;
477 inside1 = (n1 > 0.) ;
478 }
479 if (!inside1) n1 = 0. ;
480 if (inside2) {
481 n2 = nbar_ent_p2(xx2) ;
482 inside2 = (n2 > 0.) ;
483 }
484 if (!inside2) n2 = 0. ;
485 nbar1.set(l,k,j,i) = n1 ;
486 nbar2.set(l,k,j,i) = n2 ;
487 }
488 }
489 }
490 }
491 }
492
493}
494
495
496// Energy density from enthalpy
497//------------------------------
498
499Cmp Eos_bifluid::ener_ent(const Cmp& ent1, const Cmp& ent2, const Cmp& delta2,
500 int nzet, int l_min) const {
501
502 Cmp ener(ent1.get_mp()) ;
503 assert(ent1.get_etat() != ETATNONDEF) ;
504 assert(ent2.get_etat() != ETATNONDEF) ;
505 assert(delta2.get_etat() != ETATNONDEF) ;
506
507 const Map* mp = ent1.get_mp() ; // Mapping
508 assert(mp == ent2.get_mp()) ;
509 assert(mp == delta2.get_mp()) ;
510
511 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
512 ener.set_etat_zero() ;
513 return ener;
514 }
515
516 ener.allocate_all() ;
517
518 const Mg3d* mg = mp->get_mg() ; // Multi-grid
519
520 int nz = mg->get_nzone() ; // total number of domains
521
522 // resu is set to zero in the other domains :
523
524 if (l_min > 0)
525 ener.annule(0, l_min-1) ;
526
527 if (l_min + nzet < nz)
528 ener.annule(l_min + nzet, nz - 1) ;
529
530 int np0 = mp->get_mg()->get_np(0) ;
531 int nt0 = mp->get_mg()->get_nt(0) ;
532 for (int l=l_min; l<l_min+nzet; l++) {
533 assert(mp->get_mg()->get_np(l) == np0) ;
534 assert(mp->get_mg()->get_nt(l) == nt0) ;
535 }
536
537 for (int k=0; k<np0; k++) {
538 for (int j=0; j<nt0; j++) {
539 bool inside = true ;
540 bool inside1 = true ;
541 bool inside2 = true ;
542 for (int l=l_min; l< l_min + nzet; l++) {
543 for (int i=0; i<mp->get_mg()->get_nr(l); i++) {
544 double xx, xx1, xx2, n1, n2 ;
545 xx1 = ent1(l,k,j,i) ;
546 xx2 = ent2(l,k,j,i) ;
547 xx = delta2(l,k,j,i) ;
548 if (inside) {
549 inside = (!nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
550 inside1 = ((n1 > 0.)&&(xx1>0.)) ;
551 inside2 = ((n2 > 0.)&&(xx2>0.)) ;
552 }
553 if (inside) {
554 ener.set(l,k,j,i) = ener_nbar_p(n1, n2, xx) ;
555 }
556 else {
557 if (inside1) {
558 n1 = nbar_ent_p1(xx1) ;
559 inside1 = (n1 > 0.) ;
560 }
561 if (!inside1) n1 = 0. ;
562 if (inside2) {
563 n2 = nbar_ent_p2(xx2) ;
564 inside2 = (n2 > 0.) ;
565 }
566 if (!inside2) n2 = 0. ;
567 ener.set(l,k,j,i) = ener_nbar_p(n1, n2, 0.) ;
568 }
569 }
570 }
571 }
572 }
573 return ener ;
574}
575
576// Pressure from enthalpies
577//-------------------------
578
579Cmp Eos_bifluid::press_ent(const Cmp& ent1, const Cmp& ent2, const Cmp& delta2,
580 int nzet, int l_min ) const {
581
582 Cmp press(ent1.get_mp()) ;
583 assert(ent1.get_etat() != ETATNONDEF) ;
584 assert(ent2.get_etat() != ETATNONDEF) ;
585 assert(delta2.get_etat() != ETATNONDEF) ;
586
587 const Map* mp = ent1.get_mp() ; // Mapping
588 assert(mp == ent2.get_mp()) ;
589
590 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
591 press.set_etat_zero() ;
592 return press;
593 }
594 press.allocate_all() ;
595
596 const Mg3d* mg = mp->get_mg() ; // Multi-grid
597
598 int nz = mg->get_nzone() ; // total number of domains
599
600 // resu is set to zero in the other domains :
601
602 if (l_min > 0)
603 press.annule(0, l_min-1) ;
604
605 if (l_min + nzet < nz)
606 press.annule(l_min + nzet, nz - 1) ;
607
608 int np0 = mp->get_mg()->get_np(0) ;
609 int nt0 = mp->get_mg()->get_nt(0) ;
610 for (int l=l_min; l<l_min+nzet; l++) {
611 assert(mp->get_mg()->get_np(l) == np0) ;
612 assert(mp->get_mg()->get_nt(l) == nt0) ;
613 }
614
615 for (int k=0; k<np0; k++) {
616 for (int j=0; j<nt0; j++) {
617 bool inside = true ;
618 bool inside1 = true ;
619 bool inside2 = true ;
620 for (int l=l_min; l< l_min + nzet; l++) {
621 for (int i=0; i<mp->get_mg()->get_nr(l); i++) {
622 double xx, xx1, xx2, n1, n2 ;
623 xx1 = ent1(l,k,j,i) ;
624 xx2 = ent2(l,k,j,i) ;
625 xx = delta2(l,k,j,i) ;
626 if (inside) {
627 inside = (!nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
628 inside1 = ((n1 > 0.)&&(xx1>0.)) ;
629 inside2 = ((n2 > 0.)&&(xx2>0.)) ;
630 }
631 if (inside)
632 press.set(l,k,j,i) = press_nbar_p(n1, n2, xx) ;
633 else {
634 if (inside1) {
635 n1 = nbar_ent_p1(xx1) ;
636 inside1 = (n1 > 0.) ;
637 }
638 if (!inside1) n1 = 0. ;
639 if (inside2) {
640 n2 = nbar_ent_p2(xx2) ;
641 inside2 = (n2 > 0.) ;
642 }
643 if (!inside2) n2 = 0. ;
644 press.set(l,k,j,i) = press_nbar_p(n1, n2, 0. ) ;
645 }
646 }
647 }
648 }
649 }
650 return press ;
651}
652
653// Generic computational routine for get_Kxx
654//------------------------------------------
655
656void Eos_bifluid::calcule(const Cmp& nbar1, const Cmp& nbar2, const Cmp&
657 x2, int nzet, int l_min, double
658 (Eos_bifluid::*fait)(double, double, double) const,
659 Cmp& resu) const {
660
661 assert(nbar1.get_etat() != ETATNONDEF) ;
662 assert(nbar2.get_etat() != ETATNONDEF) ;
663 assert(x2.get_etat() != ETATNONDEF) ;
664
665 const Map* mp = nbar1.get_mp() ; // Mapping
666 assert(mp == nbar2.get_mp()) ;
667
668
669 if ((nbar1.get_etat() == ETATZERO)&&(nbar2.get_etat() == ETATZERO)) {
670 resu.set_etat_zero() ;
671 return ;
672 }
673
674 bool nb1 = nbar1.get_etat() == ETATQCQ ;
675 bool nb2 = nbar2.get_etat() == ETATQCQ ;
676 bool bx2 = x2.get_etat() == ETATQCQ ;
677 const Valeur* vnbar1 = 0x0 ;
678 const Valeur* vnbar2 = 0x0 ;
679 const Valeur* vx2 = 0x0 ;
680 if (nb1) { vnbar1 = &nbar1.va ;
681 vnbar1->coef_i() ; }
682 if (nb2) { vnbar2 = &nbar2.va ;
683 vnbar2->coef_i() ; }
684 if (bx2) {vx2 = & x2.va ;
685 vx2->coef_i() ; }
686
687 const Mg3d* mg = mp->get_mg() ; // Multi-grid
688
689 int nz = mg->get_nzone() ; // total number of domains
690
691 // Preparations for a point by point computation:
692 resu.set_etat_qcq() ;
693 Valeur& vresu = resu.va ;
694 vresu.set_etat_c_qcq() ;
695 vresu.c->set_etat_qcq() ;
696
697 // Loop on domains where the computation has to be done :
698 for (int l = l_min; l< l_min + nzet; l++) {
699
700 assert(l>=0) ;
701 assert(l<nz) ;
702
703 Tbl* tnbar1 = 0x0 ;
704 Tbl* tnbar2 = 0x0 ;
705 Tbl* tx2 = 0x0 ;
706
707 if (nb1) tnbar1 = vnbar1->c->t[l] ;
708 if (nb2) tnbar2 = vnbar2->c->t[l] ;
709 if (bx2) tx2 = vx2->c->t[l] ;
710 Tbl* tresu = vresu.c->t[l] ;
711
712 bool nb1b = false ;
713 if (nb1) nb1b = tnbar1->get_etat() == ETATQCQ ;
714 bool nb2b = false ;
715 if (nb2) nb2b = tnbar2->get_etat() == ETATQCQ ;
716 bool bx2b = false ;
717 if (bx2) bx2b = tx2->get_etat() == ETATQCQ ;
718 tresu->set_etat_qcq() ;
719
720 double n1, n2, xx ;
721
722 for (int i=0; i<tnbar1->get_taille(); i++) {
723
724 n1 = nb1b ? tnbar1->t[i] : 0 ;
725 n2 = nb2b ? tnbar2->t[i] : 0 ;
726 xx = bx2b ? tx2->t[i] : 0 ;
727 tresu->t[i] = (this->*fait)(n1, n2, xx ) ;
728 }
729
730
731
732 } // End of the loop on domains where the computation had to be done
733
734 // resu is set to zero in the other domains :
735
736 if (l_min > 0) {
737 resu.annule(0, l_min-1) ;
738 }
739
740 if (l_min + nzet < nz) {
741 resu.annule(l_min + nzet, nz - 1) ;
742 }
743}
744
745Cmp Eos_bifluid::get_Knn(const Cmp& nbar1, const Cmp& nbar2, const Cmp&
746 delta2, int nzet, int l_min) const {
747
748 Cmp resu(nbar1.get_mp()) ;
749
750 calcule(nbar1, nbar2, delta2, nzet, l_min, &Eos_bifluid::get_K11, resu) ;
751
752 return resu ;
753
754}
755
756Cmp Eos_bifluid::get_Knp(const Cmp& nbar1, const Cmp& nbar2, const Cmp&
757 delta2, int nzet, int l_min) const {
758
759 Cmp resu(delta2.get_mp()) ;
760
761 calcule(nbar1, nbar2, delta2, nzet, l_min, &Eos_bifluid::get_K12, resu) ;
762
763 return resu ;
764
765}
766
767Cmp Eos_bifluid::get_Kpp(const Cmp& nbar1, const Cmp& nbar2, const Cmp&
768 delta2, int nzet, int l_min) const {
769
770 Cmp resu(nbar2.get_mp()) ;
771
772 calcule(nbar1, nbar2, delta2, nzet, l_min, &Eos_bifluid::get_K22, resu) ;
773
774 return resu ;
775
776}
777
778// new functions --> ok ?
779
780// Generic computational routine for get_Kxx_ent
781//----------------------------------------------
782//Rk : nbar1 and nbar2 have been replaced by ent1 and ent2
783void Eos_bifluid::calcule_ent(const Cmp& ent1, const Cmp& ent2, const Cmp&
784 x2, int nzet, int l_min, double
785 (Eos_bifluid::*fait)(double, double, double) const,
786 Cmp& resu) const {
787
788 assert(ent1.get_etat() != ETATNONDEF) ;
789 assert(ent2.get_etat() != ETATNONDEF) ;
790 assert(x2.get_etat() != ETATNONDEF) ;
791
792 const Map* mp = ent1.get_mp() ; // Mapping
793 assert(mp == ent2.get_mp()) ;
794
795
796 if ((ent1.get_etat() == ETATZERO)&&(ent2.get_etat() == ETATZERO)) {
797 resu.set_etat_zero() ;
798 return ;
799 }
800
801 bool le1 = ent1.get_etat() == ETATQCQ ;
802 bool le2 = ent2.get_etat() == ETATQCQ ;
803 bool bx2 = x2.get_etat() == ETATQCQ ;
804 const Valeur* vent1 = 0x0 ;
805 const Valeur* vent2 = 0x0 ;
806 const Valeur* vx2 = 0x0 ;
807 if (le1) { vent1 = &ent1.va ;
808 vent1->coef_i() ; }
809 if (le2) { vent2 = &ent2.va ;
810 vent2->coef_i() ; }
811 if (bx2) {vx2 = & x2.va ;
812 vx2->coef_i() ; }
813
814 const Mg3d* mg = mp->get_mg() ; // Multi-grid
815
816 int nz = mg->get_nzone() ; // total number of domains
817
818 // Preparations for a point by point computation:
819 resu.set_etat_qcq() ;
820 Valeur& vresu = resu.va ;
821 vresu.set_etat_c_qcq() ;
822 vresu.c->set_etat_qcq() ;
823
824 // Loop on domains where the computation has to be done :
825 for (int l = l_min; l< l_min + nzet; l++) {
826
827 assert(l>=0) ;
828 assert(l<nz) ;
829
830 Tbl* tent1 = 0x0 ;
831 Tbl* tent2 = 0x0 ;
832 Tbl* tx2 = 0x0 ;
833
834 if (le1) tent1 = vent1->c->t[l] ;
835 if (le2) tent2 = vent2->c->t[l] ;
836 if (bx2) tx2 = vx2->c->t[l] ;
837 Tbl* tresu = vresu.c->t[l] ;
838
839 bool le1b = false ;
840 if (le1) le1b = tent1->get_etat() == ETATQCQ ;
841 bool le2b = false ;
842 if (le2) le2b = tent2->get_etat() == ETATQCQ ;
843 bool bx2b = false ;
844 if (bx2) bx2b = tx2->get_etat() == ETATQCQ ;
845 tresu->set_etat_qcq() ;
846
847 double H1, H2, xx ;
848
849 for (int i=0; i<tent1->get_taille(); i++) {
850
851 H1 = le1b ? tent1->t[i] : 0 ;
852 H2 = le2b ? tent2->t[i] : 0 ;
853 xx = bx2b ? tx2->t[i] : 0 ;
854
855//cout << " LA" << H1 << " " << H2 << " " << xx << endl;
856 tresu->t[i] = (this->*fait)(xx, H1, H2) ;
857 }
858
859 } // End of the loop on domains where the computation had to be done
860
861 // resu is set to zero in the other domains :
862
863 if (l_min > 0) {
864 resu.annule(0, l_min-1) ;
865 }
866
867 if (l_min + nzet < nz) {
868 resu.annule(l_min + nzet, nz - 1) ;
869 }
870}
871
872// rajout de ces fonctions
873Cmp Eos_bifluid::get_Knn_ent(const Cmp& ent1, const Cmp& ent2, const Cmp&
874 delta2, int nzet, int l_min) const {
875
876 Cmp resu(ent1.get_mp()) ;
877
878 calcule_ent(ent1, ent2, delta2, nzet, l_min, &Eos_bifluid::get_K11, resu) ;
879
880 return resu ;
881
882}
883
884Cmp Eos_bifluid::get_Knp_ent(const Cmp& ent1, const Cmp& ent2, const Cmp&
885 delta2, int nzet, int l_min) const {
886
887 Cmp resu(delta2.get_mp()) ;
888
889 calcule_ent(ent1, ent2, delta2, nzet, l_min, &Eos_bifluid::get_K12, resu) ;
890
891 return resu ;
892
893}
894
895Cmp Eos_bifluid::get_Kpp_ent(const Cmp& ent1, const Cmp& ent2, const Cmp&
896 delta2, int nzet, int l_min) const {
897
898 Cmp resu(ent2.get_mp()) ;
899
900 calcule_ent(ent1, ent2, delta2, nzet, l_min, &Eos_bifluid::get_K22, resu) ;
901
902 return resu ;
903
904}
905
906}
907
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition cmp.h:446
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition cmp.C:323
int get_etat() const
Returns the logical state.
Definition cmp.h:899
Valeur va
The numerical value of the Cmp
Definition cmp.h:464
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition cmp.C:348
Tbl & set(int l)
Read/write of the value in a given domain.
Definition cmp.h:724
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition cmp.C:289
const Map * get_mp() const
Returns the mapping.
Definition cmp.h:901
Analytic equation of state for two fluids (Newtonian case).
2-fluids equation of state base class.
void calcule(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min, double(Eos_bifluid::*fait)(double, double, double) const, Cmp &resu) const
General computational method for Cmp 's ( 's).
virtual void calcule_tout(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, Cmp &ener, Cmp &press, int nzet, int l_min=0) const
General computational method for Cmp 's, it computes both baryon densities, energy and pressure profi...
virtual double nbar_ent_p2(const double ent2) const =0
Computes baryon density out of the log-enthalpy assuming that only fluid 2 is present (virtual functi...
double get_m1() const
Return the individual particule mass
virtual void sauve(FILE *) const
Save in a file.
virtual ~Eos_bifluid()
Destructor.
void nbar_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, Cmp &nbar1, Cmp &nbar2, int nzet, int l_min=0) const
Computes both baryon density fields from the log-enthalpy fields and the relative velocity.
double m_1
Individual particle mass [unit: ].
virtual double nbar_ent_p1(const double ent1) const =0
Computes baryon density out of the log-enthalpy asuming that only fluid 1 is present (virtual functio...
string name
EOS name.
double get_m2() const
Return the individual particule mass
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to.
virtual double ener_nbar_p(const double nbar1, const double nbar2, const double delta2) const =0
Computes the total energy density from the baryonic densities and the relative velocity.
string get_name() const
Returns the EOS name.
virtual double press_nbar_p(const double nbar1, const double nbar2, const double delta2) const =0
Computes the pressure from the baryonic densities and the relative velocity.
Cmp press_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, int nzet, int l_min=0) const
Computes the pressure from the log-enthalpy fields and the relative velocity.
Cmp ener_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &delta2, int nzet, int l_min=0) const
Computes the total energy density from the log-enthalpy fields and the relative velocity.
virtual bool nbar_ent_p(const double ent1, const double ent2, const double delta2, double &nbar1, double &nbar2) const =0
Computes both baryon densities from the log-enthalpies (virtual function implemented in the derived c...
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
double m_2
Individual particle mass [unit: ].
Eos_bifluid()
Standard constructor.
Base class for coordinate mappings.
Definition map.h:670
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Definition map.h:765
Multi-domain grid.
Definition grilles.h:273
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
Definition grilles.h:462
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
Definition grilles.h:457
int get_nzone() const
Returns the number of domains.
Definition grilles.h:448
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Definition grilles.h:452
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
Definition mtbl.h:132
double * t
The array of double.
Definition tbl.h:173
Values and coefficients of a (real-value) function.
Definition valeur.h:287
Mtbl * c
Values of the function at the points of the multi-grid
Definition valeur.h:299
void coef_i() const
Computes the physical value of *this.
int fread_be(int *aa, int size, int nb, FILE *fich)
Reads integer(s) from a binary file according to the big endian convention.
Definition fread_be.C:69
int read_variable(const char *fname, const char *var_name, char *fmt, void *varp)
Reads a variable from file.
int fwrite_be(const int *aa, int size, int nb, FILE *fich)
Writes integer(s) into a binary file according to the big endian convention.
Definition fwrite_be.C:70
Lorene prototypes.
Definition app_hor.h:64