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 $" ;
136#include "eos_bifluid.h"
138#include "utilitaires.h"
149 name(
"EoS bi-fluid"), m_1(1), m_2(1)
155 name(name_i), m_1(mass1), m_2(mass2)
161 name(eos_i.name), m_1(eos_i.m_1), m_2(eos_i.m_2)
168 char dummy [MAX_EOSNAME];
169 if (fread(dummy,
sizeof(
char),MAX_EOSNAME, fich) == 0)
170 cerr <<
"Reading problem in " << __FILE__ << endl ;
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") ;
196 cerr <<
"ERROR: Eos_bifluid(const char*): could not read either of \
197the variables 'name', 'm_1, or 'm_2' from file " << fname << endl;
213 fich.getline(blabla, 80) ;
215 fich >>
m_1 ; fich.getline(blabla, 80) ;
216 fich >>
m_2 ; fich.getline(blabla, 80) ;
247 assert(
name.size() < MAX_EOSNAME) ;
248 char dummy [MAX_EOSNAME];
251 fwrite_be(&ident,
sizeof(
int), 1, fich) ;
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 ;
266 ost <<
" Mean particle 1 mass : " << eqetat.
get_m1() <<
" m_B" << endl ;
267 ost <<
" Mean particle 2 mass : " << eqetat.
get_m2() <<
" m_B" << endl ;
282 const Cmp& delta2,
Cmp& nbar1,
Cmp& nbar2,
284 int nzet,
int l_min)
const {
286 assert(ent1.
get_etat() != ETATNONDEF) ;
287 assert(ent2.
get_etat() != ETATNONDEF) ;
288 assert(delta2.
get_etat() != ETATNONDEF) ;
291 assert(mp == ent2.
get_mp()) ;
292 assert(mp == delta2.
get_mp()) ;
293 assert(mp == ener.
get_mp()) ;
314 nbar1.
annule(0, l_min-1) ;
315 nbar2.
annule(0, l_min-1) ;
317 press.
annule(0, l_min-1) ;
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) ;
329 for (
int l=l_min; l<l_min+nzet; l++) {
337 bool slow_rot_style =
false;
342 if (this_eos -> get_typeos() == 5)
344 slow_rot_style =
true;
345 cout <<
"DEBUG: using EOS-inversion slow-rot-style!! " << endl;
351 for (
int k=0; k<np0; k++) {
352 for (
int j=0; j<nt0; j++) {
354 bool inside1 = true ;
355 bool inside2 = true ;
356 for (
int l=l_min; l< l_min + nzet; l++) {
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) ;
363 inside = (!
nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
366 inside1 = (n1 > 0.) ;
368 inside2 = (n2 > 0.) ;
374 n1 = (n1 > 0) ? n1: 0;
375 n2 = (n2 > 0) ? n2: 0;
379 nbar1.
set(l,k,j,i) = n1 ;
380 nbar2.
set(l,k,j,i) = n2 ;
385 inside1 = (n1 > 0.) ;
389 inside2 = (n2 > 0.) ;
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 ;
412 Cmp& nbar1,
Cmp& nbar2,
int nzet,
int l_min)
const {
414 assert(ent1.
get_etat() != ETATNONDEF) ;
415 assert(ent2.
get_etat() != ETATNONDEF) ;
416 assert(delta2.
get_etat() != ETATNONDEF) ;
419 assert(mp == ent2.
get_mp()) ;
420 assert(mp == delta2.
get_mp()) ;
421 assert(mp == nbar1.
get_mp()) ;
438 nbar1.
annule(0, l_min-1) ;
439 nbar2.
annule(0, l_min-1) ;
442 if (l_min + nzet < nz) {
443 nbar1.
annule(l_min + nzet, nz - 1) ;
444 nbar2.
annule(l_min + nzet, nz - 1) ;
449 for (
int l=l_min; l<l_min+nzet; l++) {
454 for (
int k=0; k<np0; k++) {
455 for (
int j=0; j<nt0; j++) {
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) ;
466 inside = (!
nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
467 inside1 = ((n1 > 0.)&&(xx1>0.)) ;
468 inside2 = ((n2 > 0.)&&(xx2>0.)) ;
471 nbar1.
set(l,k,j,i) = n1 ;
472 nbar2.
set(l,k,j,i) = n2 ;
477 inside1 = (n1 > 0.) ;
479 if (!inside1) n1 = 0. ;
482 inside2 = (n2 > 0.) ;
484 if (!inside2) n2 = 0. ;
485 nbar1.
set(l,k,j,i) = n1 ;
486 nbar2.
set(l,k,j,i) = n2 ;
500 int nzet,
int l_min)
const {
503 assert(ent1.
get_etat() != ETATNONDEF) ;
504 assert(ent2.
get_etat() != ETATNONDEF) ;
505 assert(delta2.
get_etat() != ETATNONDEF) ;
508 assert(mp == ent2.
get_mp()) ;
509 assert(mp == delta2.
get_mp()) ;
527 if (l_min + nzet < nz)
528 ener.
annule(l_min + nzet, nz - 1) ;
532 for (
int l=l_min; l<l_min+nzet; l++) {
537 for (
int k=0; k<np0; k++) {
538 for (
int j=0; j<nt0; j++) {
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) ;
549 inside = (!
nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
550 inside1 = ((n1 > 0.)&&(xx1>0.)) ;
551 inside2 = ((n2 > 0.)&&(xx2>0.)) ;
559 inside1 = (n1 > 0.) ;
561 if (!inside1) n1 = 0. ;
564 inside2 = (n2 > 0.) ;
566 if (!inside2) n2 = 0. ;
580 int nzet,
int l_min )
const {
583 assert(ent1.
get_etat() != ETATNONDEF) ;
584 assert(ent2.
get_etat() != ETATNONDEF) ;
585 assert(delta2.
get_etat() != ETATNONDEF) ;
588 assert(mp == ent2.
get_mp()) ;
603 press.
annule(0, l_min-1) ;
605 if (l_min + nzet < nz)
606 press.
annule(l_min + nzet, nz - 1) ;
610 for (
int l=l_min; l<l_min+nzet; l++) {
615 for (
int k=0; k<np0; k++) {
616 for (
int j=0; j<nt0; j++) {
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) ;
627 inside = (!
nbar_ent_p(xx1, xx2, xx, n1, n2)) ;
628 inside1 = ((n1 > 0.)&&(xx1>0.)) ;
629 inside2 = ((n2 > 0.)&&(xx2>0.)) ;
636 inside1 = (n1 > 0.) ;
638 if (!inside1) n1 = 0. ;
641 inside2 = (n2 > 0.) ;
643 if (!inside2) n2 = 0. ;
657 x2,
int nzet,
int l_min,
double
658 (
Eos_bifluid::*fait)(double, double, double) const,
661 assert(nbar1.get_etat() != ETATNONDEF) ;
662 assert(nbar2.
get_etat() != ETATNONDEF) ;
663 assert(x2.
get_etat() != ETATNONDEF) ;
666 assert(mp == nbar2.
get_mp()) ;
670 resu.set_etat_zero() ;
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 ;
680 if (nb1) { vnbar1 = &nbar1.
va ;
682 if (nb2) { vnbar2 = &nbar2.
va ;
684 if (bx2) {vx2 = & x2.
va ;
687 const Mg3d* mg = mp->
get_mg() ;
692 resu.set_etat_qcq() ;
693 Valeur& vresu = resu.va ;
694 vresu.set_etat_c_qcq() ;
695 vresu.c->set_etat_qcq() ;
698 for (
int l = l_min; l< l_min + nzet; l++) {
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] ;
713 if (nb1) nb1b = tnbar1->get_etat() == ETATQCQ ;
715 if (nb2) nb2b = tnbar2->get_etat() == ETATQCQ ;
717 if (bx2) bx2b = tx2->get_etat() == ETATQCQ ;
718 tresu->set_etat_qcq() ;
722 for (
int i=0; i<tnbar1->get_taille(); i++) {
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 ) ;
737 resu.annule(0, l_min-1) ;
740 if (l_min + nzet < nz) {
741 resu.annule(l_min + nzet, nz - 1) ;
746 delta2,
int nzet,
int l_min)
const {
750 calcule(nbar1, nbar2, delta2, nzet, l_min, &Eos_bifluid::get_K11, resu) ;
757 delta2,
int nzet,
int l_min)
const {
761 calcule(nbar1, nbar2, delta2, nzet, l_min, &Eos_bifluid::get_K12, resu) ;
768 delta2,
int nzet,
int l_min)
const {
772 calcule(nbar1, nbar2, delta2, nzet, l_min, &Eos_bifluid::get_K22, resu) ;
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,
788 assert(ent1.get_etat() != ETATNONDEF) ;
789 assert(ent2.
get_etat() != ETATNONDEF) ;
790 assert(x2.
get_etat() != ETATNONDEF) ;
793 assert(mp == ent2.
get_mp()) ;
797 resu.set_etat_zero() ;
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 ;
807 if (le1) { vent1 = &ent1.
va ;
809 if (le2) { vent2 = &ent2.
va ;
811 if (bx2) {vx2 = & x2.
va ;
814 const Mg3d* mg = mp->
get_mg() ;
819 resu.set_etat_qcq() ;
820 Valeur& vresu = resu.va ;
821 vresu.set_etat_c_qcq() ;
822 vresu.c->set_etat_qcq() ;
825 for (
int l = l_min; l< l_min + nzet; l++) {
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] ;
840 if (le1) le1b = tent1->get_etat() == ETATQCQ ;
842 if (le2) le2b = tent2->get_etat() == ETATQCQ ;
844 if (bx2) bx2b = tx2->get_etat() == ETATQCQ ;
845 tresu->set_etat_qcq() ;
849 for (
int i=0; i<tent1->get_taille(); i++) {
851 H1 = le1b ? tent1->t[i] : 0 ;
852 H2 = le2b ? tent2->t[i] : 0 ;
853 xx = bx2b ? tx2->t[i] : 0 ;
856 tresu->t[i] = (this->*fait)(xx, H1, H2) ;
864 resu.annule(0, l_min-1) ;
867 if (l_min + nzet < nz) {
868 resu.annule(l_min + nzet, nz - 1) ;
873Cmp Eos_bifluid::get_Knn_ent(
const Cmp& ent1,
const Cmp& ent2,
const Cmp&
874 delta2,
int nzet,
int l_min)
const {
878 calcule_ent(ent1, ent2, delta2, nzet, l_min, &Eos_bifluid::get_K11, resu) ;
884Cmp Eos_bifluid::get_Knp_ent(
const Cmp& ent1,
const Cmp& ent2,
const Cmp&
885 delta2,
int nzet,
int l_min)
const {
889 calcule_ent(ent1, ent2, delta2, nzet, l_min, &Eos_bifluid::get_K12, resu) ;
895Cmp Eos_bifluid::get_Kpp_ent(
const Cmp& ent1,
const Cmp& ent2,
const Cmp&
896 delta2,
int nzet,
int l_min)
const {
900 calcule_ent(ent1, ent2, delta2, nzet, l_min, &Eos_bifluid::get_K22, resu) ;
Component of a tensorial field *** DEPRECATED : use class Scalar instead ***.
void allocate_all()
Sets the logical state to ETATQCQ (ordinary state) and performs the memory allocation of all the elem...
int get_etat() const
Returns the logical state.
Valeur va
The numerical value of the Cmp
void annule(int l)
Sets the Cmp to zero in a given domain.
Tbl & set(int l)
Read/write of the value in a given domain.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
const Map * get_mp() const
Returns the mapping.
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...
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.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's.
double * t
The array of double.
Values and coefficients of a (real-value) function.
Mtbl * c
Values of the function at the points of the multi-grid
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.
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.