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 
29 char 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 // ---------------------------------
147 namespace Lorene {
149  name("EoS bi-fluid"), m_1(1), m_2(1)
150 { }
151 
152 // Standard constructor with name and masses
153 // ---------------------------------
154 Eos_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 // ---------------------------------
179 Eos_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 \
197 the 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 // ---------------------------------
205  Eos_bifluid::Eos_bifluid(ifstream& fich){
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 
227 { }
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 
245 void 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 
264 ostream& 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 
281 void 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 
411 void 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 
499 Cmp 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 
579 Cmp 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 
656 void 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 
745 Cmp 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 
756 Cmp 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 
767 Cmp 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
783 void 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
873 Cmp 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 
884 Cmp 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 
895 Cmp 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 
Cmp get_Kpp(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 2) .
Definition: eos_bifluid.C:767
virtual double get_K12(const double n1, const double n2, const double x) const =0
Computes the derivative of the energy with respect to .
const Map * get_mp() const
Returns the mapping.
Definition: cmp.h:901
double get_m2() const
Return the individual particule mass .
Definition: eos_bifluid.h:262
virtual void sauve(FILE *) const
Save in a file.
Definition: eos_bifluid.C:245
virtual double get_K22(const double n1, const double n2, const double x) const =0
Computes the derivative of the energy/(baryonic density 2) .
double m_1
Individual particle mass [unit: ].
Definition: eos_bifluid.h:184
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
Definition: cmp.h:446
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.
Definition: eos_bifluid.C:579
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...
void annule(int l)
Sets the Cmp to zero in a given domain.
Definition: cmp.C:348
Cmp get_Knn_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 1) .
Definition: eos_bifluid.C:873
Lorene prototypes.
Definition: app_hor.h:64
string name
EOS name.
Definition: eos_bifluid.h:179
virtual int identify() const =0
Returns a number to identify the sub-classe of Eos_bifluid the object belongs to. ...
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.
Definition: eos_bifluid.C:499
int get_etat() const
Returns the logical state.
Definition: cmp.h:899
void coef_i() const
Computes the physical value of *this.
Base class for coordinate mappings.
Definition: map.h:670
Values and coefficients of a (real-value) function.
Definition: valeur.h:287
int get_etat() const
Gives the logical state.
Definition: tbl.h:394
void operator=(const Eos_bifluid &)
Assignment to another Eos_bifluid.
Definition: eos_bifluid.C:232
2-fluids equation of state base class.
Definition: eos_bifluid.h:173
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: tbl.C:361
Cmp get_Knn(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 1) .
Definition: eos_bifluid.C:745
virtual double get_K11(const double n1, const double n2, const double x) const =0
Computes the derivative of the energy with respect to (baryonic density 1) .
double m_2
Individual particle mass [unit: ].
Definition: eos_bifluid.h:189
Cmp get_Kpp_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy/(baryonic density 2) .
Definition: eos_bifluid.C:895
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
Definition: cmp.C:289
Eos_bifluid()
Standard constructor.
Definition: eos_bifluid.C:148
Cmp get_Knp_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy with respect to .
Definition: eos_bifluid.C:884
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.
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...
Definition: eos_bifluid.C:411
double * t
The array of double.
Definition: tbl.h:173
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...
Mtbl * c
Values of the function at the points of the multi-grid.
Definition: valeur.h:299
int get_nzone() const
Returns the number of domains.
Definition: grilles.h:448
Analytic equation of state for two fluids (Newtonian case).
Definition: eos_bifluid.h:1258
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
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: cmp.C:304
friend ostream & operator<<(ostream &, const Eos_bifluid &)
Display.
Definition: eos_bifluid.C:264
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition: mtbl.C:299
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...
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 &#39;s ( &#39;s).
Definition: eos_bifluid.C:656
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.
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
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
Definition: cmp.C:323
Multi-domain grid.
Definition: grilles.h:273
Tbl & set(int l)
Read/write of the value in a given domain.
Definition: cmp.h:724
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 &#39;s, it computes both baryon densities, energy and pressure profi...
Definition: eos_bifluid.C:281
virtual ~Eos_bifluid()
Destructor.
Definition: eos_bifluid.C:226
string get_name() const
Returns the EOS name.
Definition: eos_bifluid.h:246
int read_variable(const char *fname, const char *var_name, char *fmt, void *varp)
Reads a variable from file.
int get_taille() const
Gives the total size (ie dim.taille)
Definition: tbl.h:397
void set_etat_c_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl c )...
Definition: valeur.C:701
Basic array class.
Definition: tbl.h:161
Cmp get_Knp(const Cmp &nbar1, const Cmp &nbar2, const Cmp &x2, int nzet, int l_min=0) const
Computes the derivatives of the energy with respect to .
Definition: eos_bifluid.C:756
Tbl ** t
Array (size nzone ) of pointers on the Tbl &#39;s.
Definition: mtbl.h:132
Valeur va
The numerical value of the Cmp.
Definition: cmp.h:464
double get_m1() const
Return the individual particule mass .
Definition: eos_bifluid.h:256
void calcule_ent(const Cmp &ent1, const Cmp &ent2, const Cmp &x2, int nzet, int l_min, double(Eos_bifluid::*fait)(double, double, double) const, Cmp &resu) const
General computational method for Cmp &#39;s ( &#39;s).
Definition: eos_bifluid.C:783