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