Commit 12fb511b authored by Francesco's avatar Francesco
Browse files

Updates2 GRMHDb (trying to resume the 6month old code)

parent 70dffd3f

ApplicationExamples/GRMHDb/GRMHDb_ADERDG_FV_2D_clean/GRMHDbSolver_ADERDG.cpp

+76 −148
Original line number Diff line number Diff line
@@ -180,6 +180,7 @@ void GRMHDb::GRMHDbSolver_ADERDG::boundaryValues(const double* const x, const do 
  stateOut[16] = 0.0;

  stateOut[17] = 0.0;

  stateOut[18] = 0.0;



  fluxOut[0] = 0.0;

  fluxOut[1] = 0.0;

  fluxOut[2] = 0.0;
@@ -200,6 +201,40 @@ void GRMHDb::GRMHDbSolver_ADERDG::boundaryValues(const double* const x, const do 
  fluxOut[17] = 0.0;

  fluxOut[18] = 0.0;





 /* stateOut[0] = exp(-(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / 8.0);

  stateOut[1] = sin(x[1])*sin(x[0]);

  stateOut[2] = sin(x[2]);

  for (int i = 3; i < numberOfData; i++) {

	  stateOut[i] = cos(x[0]);

  }

  fluxOut[0] = exp(-(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / 8.0);

  fluxOut[1] = sin(x[1])*sin(x[0]);

  fluxOut[2] = sin(x[2]);

  for (int i = 3; i < numberOfData; i++) {

	  fluxOut[i] = cos(x[0]);

  }*/

 

 

 //// STANDARD ANALYTICAL BOUNDARY CONDITIONS:



 for(int dd=0; dd<nDim; dd++) F[dd] = Fs[dd];



 for(int i=0; i < basisSize; i++)  { // i == time

	  const double weight = kernels::legendre::weights[order][i];

	  const double xi = kernels::legendre::nodes[order][i];

	  double ti = t + xi * dt;



	  initialdata_(x, &ti, Qgp);

	  //pdeflux_(F[0], F[1], F[2], Qgp);

	  flux(Qgp, F);

          for(int m=0; m < numberOfData; m++) {

		  stateOut[m] += weight * Qgp[m];

		  fluxOut[m] += weight * Fs[normalNonZero][m];

	  }

 }



/*

  // THIS IS FOR THE 1D Riemann problems (invisicd reflective boundary conditions at the y boundaries)

  if (normalNonZero == 0) {

	  // STANDARD ANALYTICAL BOUNDARY CONDITIONS:
@@ -225,46 +260,7 @@ void GRMHDb::GRMHDbSolver_ADERDG::boundaryValues(const double* const x, const do 
	  flux(stateOut, F); 

	  std::copy_n(F[normalNonZero], numberOfVariables, fluxOut);

  }



  //printf("\n I WAS HERE!!!!!!!!");



 //// STANDARD ANALYTICAL BOUNDARY CONDITIONS:

 //for(int dd=0; dd<nDim; dd++) F[dd] = Fs[dd];



 //for(int i=0; i < basisSize; i++)  { // i == time

	//  const double weight = kernels::legendre::weights[order][i];

	//  const double xi = kernels::legendre::nodes[order][i];

	//  double ti = t + xi * dt;



	//  initialdata_(x, &ti, Qgp);

	//  //pdeflux_(F[0], F[1], F[2], Qgp);

	//  flux(Qgp, F);

 //         for(int m=0; m < numberOfData; m++) {

	//	  stateOut[m] += weight * Qgp[m];

	//	  fluxOut[m] += weight * Fs[normalNonZero][m];

	//  }

 //}









 /* stateOut[0] = exp(-(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / 8.0);

  stateOut[1] = sin(x[1])*sin(x[0]);

  stateOut[2] = sin(x[2]);

  for (int i = 3; i < numberOfData; i++) {

	  stateOut[i] = cos(x[0]);

  }

  fluxOut[0] = exp(-(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / 8.0);

  fluxOut[1] = sin(x[1])*sin(x[0]);

  fluxOut[2] = sin(x[2]);

  for (int i = 3; i < numberOfData; i++) {

	  fluxOut[i] = cos(x[0]);

  }*/

  /*

   double x3D[3] = {0.};

   for(int i=0;i<DIMENSIONS;i++){

	 x3D[i]=x[i];

	 }*/

*/



  /*

	for(int m=0; m < numberOfData; m++) {
@@ -272,9 +268,6 @@ void GRMHDb::GRMHDbSolver_ADERDG::boundaryValues(const double* const x, const do 
	fluxOut[m] = fluxIn[m];

	}

	*/

	

       // printf("\n******* DONE*****************");



}



exahype::solvers::Solver::RefinementControl GRMHDb::GRMHDbSolver_ADERDG::refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& center,const tarch::la::Vector<DIMENSIONS,double>& dx,double t,const int level) {
@@ -441,9 +434,6 @@ void GRMHDb::GRMHDbSolver_ADERDG::flux(const double* const Q,double** const F) { 
  // Dimensions                        = 3

  // Number of variables + parameters  = 19 + 0

  

       // printf("\n*******   *****************");



	//printf("\nSONO QUI IN flux ADERDG");

  // @todo Please implement/augment if required

  F[0][0] = 0.0;

  F[0][1] = 0.0;
@@ -520,39 +510,6 @@ void GRMHDb::GRMHDbSolver_ADERDG::flux(const double* const Q,double** const F) { 


  //pdeflux_(&F[0][0], Q);

 



 // F[1][0] = 0.0;

 // F[1][1] = 0.0;

 // F[1][2] = 0.0;

 // F[1][3] = 0.0;

 // F[1][4] = 0.0;

 // F[1][5] = 0.0;

 // F[1][6] = 0.0;

 // F[1][7] = 0.0;

 // F[1][8] = 0.0;

 // F[1][9] = 0.0;

 // F[1][10] = 0.0;

 // F[1][11] = 0.0;

 // F[1][12] = 0.0;

 // F[1][13] = 0.0;

 // F[1][14] = 0.0;

 // F[1][15] = 0.0;

 // F[1][16] = 0.0;

 // F[1][17] = 0.0;

 // F[1][18] = 0.0;





 // printf("\nSONO QUI IN FLUX ADERDG");

 ////   if(DIMENSIONS == 2){

	////	constexpr int numberOfVariables = AbstractGRMHDbSolver_ADERDG::NumberOfVariables;

	////	constexpr int numberOfParameters = AbstractGRMHDbSolver_ADERDG::NumberOfParameters;

	////	constexpr int numberOfData = numberOfVariables + numberOfParameters;

	////	double F_3[numberOfData];

	////	pdeflux_(F[0], F[1],F_3, Q);

	////}else{

	////	pdeflux_(F[0], F[1],F[2], Q);

	////}

   

}
@@ -562,7 +519,6 @@ void GRMHDb::GRMHDbSolver_ADERDG::nonConservativeProduct(const double* const Q, 
  // Tip: See header file "GRMHDb::AbstractGRMHDbSolver_ADERDG.h" for toolkit generated compile-time 

  //      constants such as Order, NumberOfVariables, and NumberOfParameters.

  

  //printf("\n ************* nonConservativeProduct ADERDG ***************");

  // @todo Please implement/augment if required

  BgradQ[0] = 0.0;

  BgradQ[1] = 0.0;
@@ -583,28 +539,7 @@ void GRMHDb::GRMHDbSolver_ADERDG::nonConservativeProduct(const double* const Q, 
  BgradQ[16] = 0.0;

  BgradQ[17] = 0.0;

  BgradQ[18] = 0.0;



  //return;

  

 /* const int numberOfVariables = AbstractGRMHDbSolver_ADERDG::NumberOfVariables;

  constexpr int tot3Dvar = 3*numberOfVariables;

  double gradQ3D[tot3Dvar]={0.};

  int count=0;

  for(int i=0;i<DIMENSIONS;i++){

   for(int m=0;m<numberOfVariables;m++){

     gradQ3D[count]=gradQ[count];

     count++;

     }

    }

    */



  pdencp_(BgradQ, Q, gradQ);





  //printf("\n *************t ADERDG ***************");







}



void GRMHDb::GRMHDbSolver_ADERDG::mapDiscreteMaximumPrincipleObservables(
@@ -781,8 +716,6 @@ void GRMHDb::GRMHDbSolver_ADERDG::riemannSolver(double* const FL, double* const 


#else







#include "kernels/aderdg/generic/Kernels.h" 

void GRMHDb::GRMHDbSolver_ADERDG::riemannSolver(double* const FL, double* const FR, const double* const QL, const double* const QR, const double t, const double dt, const tarch::la::Vector<DIMENSIONS, double>& dx, const int direction, bool isBoundaryFace, int faceIndex) {

	assertion2(direction >= 0, dt, direction);
@@ -827,10 +760,5 @@ void GRMHDb::GRMHDbSolver_ADERDG::riemannSolver(double* const FL, double* const 




}









#endif

ApplicationExamples/GRMHDb/GRMHDb_ADERDG_FV_2D_clean/GRMHDbSolver_FV.cpp

+20 −140
Original line number Diff line number Diff line
@@ -84,7 +84,6 @@ void GRMHDb::GRMHDbSolver_FV::adjustSolution(const double* const x,const double 
	// @todo Please implement/augment if required

  if (tarch::la::equals(t,0.0)) {

	  //const int nVar = GRMHDb::AbstractGRMHDbSolver_ADERDG::NumberOfVariables;

	  

	  constexpr int numberOfVariables = AbstractGRMHDbSolver_FV::NumberOfVariables;

	  constexpr int numberOfParameters = AbstractGRMHDbSolver_FV::NumberOfParameters;

	  constexpr int numberOfData = numberOfVariables + numberOfParameters;
@@ -114,15 +113,8 @@ void GRMHDb::GRMHDbSolver_FV::adjustSolution(const double* const x,const double 
		// everything in here is thread-safe w.r.t. the lock

		// call Fortran routines

		/***********************/

/*

        double x3D[3]={0.0};

        for(int i=0;i<DIMENSIONS;i++){

        x3D[i]=x[i];

        }

        //printf("x3d_FV:  %f,  %f",x3D[0],x3D[1]);

        */



        initialdata_(x, &t, Q);

       // printf("\nx,   rho:  %f,  %f",x3D[0],Q[1]);



		/************/

		lock.free();
@@ -148,12 +140,13 @@ void GRMHDb::GRMHDbSolver_FV::eigenvalues(const double* const Q, const int dInde 
  // Number of variables    = 19 + #parameters



  //printf("\n******* EIGENVALUES FV *****************");

  /*

	constexpr int numberOfVariables = AbstractGRMHDbSolver_FV::NumberOfVariables;

	constexpr int numberOfParameters = AbstractGRMHDbSolver_FV::NumberOfParameters;

	constexpr int numberOfData = numberOfVariables + numberOfParameters;

	for (int i = 0; i < numberOfData; i++) {

		assertion(!std::isnan(Q[i]));

     }

     }*/

  // @todo Please implement/augment if required

  lambda[0]  = 0.0;

  lambda[1]  = 0.0;
@@ -221,8 +214,14 @@ void GRMHDb::GRMHDbSolver_FV::boundaryValues( 
  stateOutside[17] = stateInside[17];

  stateOutside[18] = stateInside[18];



  //return;

  // d is one between 0,1,2



// THIS IS FOR ANALYTICAL BOUNDARY CONDITIONS:

  double ti = t + 0.5 * dt;

  initialdata_(x, &ti, &Qgp[0]);

  for(int m=0; m < numberOfData; m++) {

        stateOutside[m] = Qgp[m];

  }

/*

  // THIS IS FOR 1D Riemann problems. (inviscid reflection at the y boundaries)

  if (d==0) {

	  double ti = t + 0.5 * dt;
@@ -246,19 +245,7 @@ void GRMHDb::GRMHDbSolver_FV::boundaryValues( 
	  pdeprim2cons_(&stateOutside[0], &VtmpOut[0]);

	  //stateOutside[1 + d] = -stateInside[1 + d];

  }



// THIS IS FOR ANALYTICAL BOUNDARY CONDITIONS:

 // double ti = t + 0.5 * dt;

 ///* 

 //       double x3D[3]={0.0};

 //       for(int i=0;i<DIMENSIONS;i++){

 //       x3D[i]=x[i];

 //       }*/



 // initialdata_(x, &ti, &Qgp[0]);

 // for(int m=0; m < numberOfData; m++) {

 //       stateOutside[m] = Qgp[m];

 // }

*/





  /*stateOutside[0] = exp(-(x[0] * x[0] + x[1] * x[1] + x[2] * x[2]) / 8.0);
@@ -288,9 +275,6 @@ void GRMHDb::GRMHDbSolver_FV::flux(const double* const Q,double** const F) { 
	constexpr int numberOfParameters = AbstractGRMHDbSolver_FV::NumberOfParameters;

	constexpr int numberOfData = numberOfVariables + numberOfParameters;

  

  

       // printf("\n******* FLUXES *****************");



  // @todo Please implement/augment if required

  F[0][0] = 0.0;

  F[0][1] = 0.0;
@@ -364,60 +348,6 @@ void GRMHDb::GRMHDbSolver_FV::flux(const double* const Q,double** const F) { 
		pdeflux_(F[0], F[1],F[2], Q);

	}

  



  //constexpr int nFF = numberOfVariables * DIMENSIONS;

  //double FF[nFF];

  //int count;

  ////count = 0;

  ////for (int dloc = 0; dloc < DIMENSIONS; dloc++) {

	 //// for (int i = 0; i < numberOfVariables; i++) {

		////  FF[count]= F[dloc][i];

		////  printf("FF(count) is  %d,    %f", count, FF[count]);

		////  count++;

	 //// }

  ////}

  //pdeflux_(&FF[0], Q);





  //count = 0;

  //for (int dloc = 0; dloc < DIMENSIONS; dloc++) {

  //for (int i = 0; i < numberOfVariables; i++) {

	 // F[dloc][i] = FF[count];

	 // count++;

	 // }

  //}



  //F[1][0] = 0.0;

  //F[1][1] = 0.0;

  //F[1][2] = 0.0;

  //F[1][3] = 0.0;

  //F[1][4] = 0.0;

  //F[1][5] = 0.0;

  //F[1][6] = 0.0;

  //F[1][7] = 0.0;

  //F[1][8] = 0.0;

  //F[1][9] = 0.0;

  //F[1][10] = 0.0;

  //F[1][11] = 0.0;

  //F[1][12] = 0.0;

  //F[1][13] = 0.0;

  //F[1][14] = 0.0;

  //F[1][15] = 0.0;

  //F[1][16] = 0.0;

  //F[1][17] = 0.0;

  //F[1][18] = 0.0;









 //   if(DIMENSIONS == 2){

 //       //const int nVar = GRMHDb::AbstractGRMHDbSolver_FV::NumberOfVariables;

	//	double F_3[numberOfData];

	//	pdeflux_(F[0], F[1],F_3, Q);

	//}else{

	//	pdeflux_(F[0], F[1],F[2], Q);

	//}

  

}
@@ -428,8 +358,6 @@ void GRMHDb::GRMHDbSolver_FV::nonConservativeProduct(const double* const Q,cons 
  // Tip: See header file "GRMHDb::AbstractGRMHDbSolver_FV.h" for toolkit generated compile-time 

  //      constants such as PatchSize, NumberOfVariables, and NumberOfParameters.

  

       // printf("\n*******nonConservativeProduct *****************");

  

  // @todo Please implement/augment if required

  BgradQ[0] = 0.0;

  BgradQ[1] = 0.0;
@@ -450,38 +378,7 @@ void GRMHDb::GRMHDbSolver_FV::nonConservativeProduct(const double* const Q,cons 
  BgradQ[16] = 0.0;

  BgradQ[17] = 0.0;

  BgradQ[18] = 0.0;

  

  //return;

  //

  ////printf("\n QUI QUO QUA");

  //const int numberOfVariables = AbstractGRMHDbSolver_ADERDG::NumberOfVariables;

  //constexpr int tot3Dvar = 3*numberOfVariables;

  ////printf("\nTOT 3D var = %d",tot3Dvar);

  //double gradQ3D[tot3Dvar]={0.0};

  //int count;

  //count = 0;

  //for(int i=0;i<DIMENSIONS;i++){

  //for(int m=0;m<numberOfVariables;m++){

  //      //printf("\n NCP count=%d    inizio",count);

  //      //fflush(stdout);

  //      gradQ3D[count]=1.0;

  //      //printf("\n NCP count=%d  mezzo",count); 

  //      gradQ3D[count]=1.0;

  //      //fflush(stdout);



  //      double myself = gradQ[count];      

  //      //printf("\n NCP count=%d  fine ",count);

  //      //fflush(stdout);

  //      gradQ3D[count]=gradQ[count];

  //      count++;

  //  }

  //  }

   //printf("\n ARRIVATO QUI");   

  pdencp_(BgradQ, Q, gradQ);





      //  printf("\n******* FV *****************");



}
@@ -489,20 +386,10 @@ void GRMHDb::GRMHDbSolver_FV::referenceSolution(const double* const x, double t, 
	constexpr int numberOfVariables = AbstractGRMHDbSolver_FV::NumberOfVariables;

	constexpr int numberOfParameters = AbstractGRMHDbSolver_FV::NumberOfParameters;

	constexpr int numberOfData = numberOfVariables + numberOfParameters;

	

	

       // printf("\n*******referenceSolution*****************");



	

	int iErr;

	double Qcons[numberOfData];

	iErr = 0;

        /*

        double x3D[3]={0.0};

        for(int i=0;i<DIMENSIONS;i++){

        x3D[i]=x[i];

        }*/

	//printf("\nSONO QUI IN REFERENCE SOLUTION");



	initialdata_(x, &t, &Qcons[0]);



	//// test:
@@ -520,12 +407,9 @@ void GRMHDb::GRMHDbSolver_FV::referenceSolution(const double* const x, double t, 
double GRMHDb::GRMHDbSolver_FV::riemannSolver(double* fL, double *fR, const double* qL, const double* qR, const double* gradQL, const double* gradQR, const double* cellSize, int direction) {

//double GRMHDb::GRMHDbSolver_FV::riemannSolver(double* const fL, double* const fR, const double* const qL, const double* const qR, const double* gradQL, const double* gradQR, int direction) {

	//// Default FV Riemann Solver

        // printf("\n*******riemannSolver(*****************");



    //return kernels::finitevolumes::riemannsolvers::c::rusanov<true, true, false, GRMHDbSolver_FV>(*static_cast<GRMHDbSolver_FV*>(this), fL,fR,qL,qR,gradQL, gradQR, cellSize, direction);

	constexpr int numberOfVariables = AbstractGRMHDbSolver_FV::NumberOfVariables;



	//printf("SONO QUI IN riemannSolver");

	/* HLLEM */

	

	//const int numberOfVariables = GRMHDb::AbstractGRMHDbSolver_FV::NumberOfVariables;
@@ -539,8 +423,6 @@ double GRMHDb::GRMHDbSolver_FV::riemannSolver(double* fL, double *fR, const doub 


										 // std::cout << "opened ---------------------"<< std::endl;



        // printf("\n******* RIEMANN SOLVER FV*****************");



	kernels::idx2 idx_QLR(basisSize, numberOfData);

	for (int j = 0; j < basisSize; j++) {

		const double weight = kernels::legendre::weights[order][j];
@@ -551,8 +433,6 @@ double GRMHDb::GRMHDbSolver_FV::riemannSolver(double* fL, double *fR, const doub 
		}

	}

	 

        // printf("\n***DONE*****");



	double lambda = 2.0;

	hllemfluxfv_(fL, fR, qL, qR, QavL, QavR, &direction);

ApplicationExamples/GRMHDb/GRMHDb_ADERDG_FV_2D_clean/InitialData.f90

+2 −4
Original line number Diff line number Diff line 
! GRMHDb Initial Data

 

#define GRMHD

#define NoRNSTOV



RECURSIVE SUBROUTINE PDESetup(myrank)

	USE, INTRINSIC :: ISO_C_BINDING

ApplicationExamples/GRMHDb/GRMHDb_ADERDG_FV_2D_clean/PDE.f90

+36 −77
Original line number Diff line number Diff line 
! GRMHD PDE.f90

! Trento (EQNTYPE4)

!!!#define DEBUGGONE





RECURSIVE SUBROUTINE PDEPrim2Cons(Q,V)

  USE Parameters, ONLY: nVar
@@ -64,37 +64,13 @@ RECURSIVE SUBROUTINE PDEFlux(f,g,h,Q) 
  !INTENT(OUT) :: FF

  ! Local varialbes

  INTEGER :: i

  REAL :: FF(nVar,d), V(nVar)

  REAL :: FF(nVar,ndim), V(nVar)

  !REAL :: V(nVar)

  REAL, PARAMETER :: epsilon = 1e-14 

  !

#ifdef DEBUGGONE

  f = 0.0

  g = 0.0

  h = 0.0

  !FF = 0.

  STOP

  RETURN

#endif

  !

!IF(ANY(Q(6:8).NE.0)) THEN

!    PRINT *, "PDEFlux Q(6:8)",Q(6:8)

!    continue

!ENDIF

!

  !PRINT *, '*************************************'

  !PRINT *, 'PDEFluxGRMHD'

  !PRINT *, FF(:,1)

  !PRINT *, FF(:,2)

  !PRINT *, nDim

  !PRINT *, '*************************************'

  CALL PDEFluxGRMHD(FF,Q)

  !FF = 0.

!IF(ANY(FF(6:8,:).NE.0)) THEN

!    PRINT *, "PDEFlux FF(6:8,1)",FF(6:8,1)

!    PRINT *, "PDEFlux FF(6:8,2)",FF(6:8,2)

!    continue

!ENDIF

  !

#if defined(Dim2)

  DO i=1,nVar
@@ -135,6 +111,12 @@ RECURSIVE SUBROUTINE PDEFlux(f,g,h,Q) 
!		PRINT *,' PDEFlux NAN ' 

!		STOP

!	ENDIF 

!#if defined(Dim3)

!	IF( ANY( ISNAN(h) )) THEN

!		PRINT *,' PDEFlux NAN ' 

!		STOP

!	ENDIF

!#endif

!	IF(ANY(Q(6:8).NE.0.0)) THEN

!        PRINT *, Q(6:8)

!        PRINT *, "I feel magnetized :-("
@@ -161,11 +143,6 @@ RECURSIVE SUBROUTINE PDENCP(BgradQ,Q,gradQ) 
   RETURN

#endif

	!

!IF(ANY(Q(6:8).NE.0)) THEN

!    PRINT *, "PDENCP Q(6:8)",Q(6:8) 

!    continue

!ENDIF

  !PRINT *, 'PDENCPGRMHD'

   CALL PDENCPGRMHD(BgradQ,Q,gradQ)



!IF(ANY(BgradQ(6:8).NE.0)) THEN
@@ -202,6 +179,12 @@ RECURSIVE SUBROUTINE PDENCP(BgradQ,Q,gradQ) 
!		PRINT *,'---------------' 

!		STOP

!	ENDIF 

!#if defined(Dim3)

!	IF( ANY( ISNAN(gradQ(:,3)) )) THEN

!		PRINT *,'gradQ, PDENCP NAN ' 

!		STOP

!	ENDIF 

!#endif

!	IF( ANY( ISNAN(BgradQ) )) THEN

!		PRINT *,'BgradQ,  PDENCP NAN ' 

!		WRITE(*,*) Q(1:5)
@@ -249,15 +232,11 @@ RECURSIVE SUBROUTINE PDEEigenvalues(L,Q,n) 
   USE GRMHD_Mod

  IMPLICIT NONE

  REAL :: L(nVar), n(nDim), Q(nVar), Vp(nVar)

  INTEGER :: i,iErr

  INTENT(IN)  :: Q,n

  INTENT(OUT) :: L 

  ! Local variables

  

!IF(ANY(Q(6:8).NE.0)) THEN

!    PRINT *, "PDEEigenvalues Q(6:8)",Q(6:8) 

!    continue

!ENDIF



  !

  L(:) = 0

  !
@@ -269,12 +248,7 @@ RECURSIVE SUBROUTINE PDEEigenvalues(L,Q,n) 
  !

  ! L = 1.0

  CALL PDEEigenvaluesGRMHD(L,Q,n)

  !CALL PDECons2PrimGRMHD(Vp,Q,iErr)

  !

  !    WRITE(*,'(a,f18.10,f18.10)') "n(1),n(2):",n(1),n(2)

  !DO i = 1, nVar

  !    WRITE(*,'(a,i9,E16.6,E16.6,E16.6)') "i,Q(i),Vp(i),L(i):",i,Q(i),Vp(i),L(i)

  !ENDDO 

	!IF( ANY( ISNAN(L) )) THEN

	!	PRINT *,' PDEEigenvalues NAN ' 

	!	STOP
@@ -365,10 +339,6 @@ RECURSIVE SUBROUTINE PDEMatrixB(An,Q,nv) 
  ! we use this only for the Roe Matrix

	! --------------------------------------------

	!

!IF(ANY(Q(6:8).NE.0)) THEN

!    PRINT *, "PDEMatrixB Q(6:8)",Q(6:8) 

!    continue

!ENDIF

	CALL PDEMatrixBGRMHD(An,Q,nv)  

	!

	continue
@@ -478,7 +448,7 @@ END SUBROUTINE RoeMatrix 
RECURSIVE SUBROUTINE HLLEMFluxFV(FL,FR,QL,QR,QavL,QavR,NormalNonZero) 

  USE Parameters, ONLY : nVar, nDim, nLin

  USE iso_c_binding 

  USE GRMHD_Mod

  IMPLICIT NONE

  ! Local variables

  INTEGER, INTENT(IN)   :: NormalNonZero

  REAL, INTENT(IN)     :: QL(nVar)
@@ -487,21 +457,21 @@ RECURSIVE SUBROUTINE HLLEMFluxFV(FL,FR,QL,QR,QavL,QavR,NormalNonZero) 
  REAL, INTENT(INOUT)  :: FR(nVar)

  REAL    :: QavL(nVar), QavR(nVar)  

    ! Local variables 

  INTEGER           :: i,j,k,l, ml(1)  ,iErr

  REAL              :: smax, Qav(nVar)

  INTEGER           :: i,j,k,l, ml(1)  

  REAL              :: smax, Qav(nVar),sL,sR

  REAL              ::  nv(nDim), flattener(nLin)

  REAL    :: absA(nVar,nVar), amax  

  REAL    :: absA(nVar,nVar), amax, minL, maxR, minM, maxM   

  REAL    :: QM(nVar),LL(nVar),LR(nVar),LM(nVar)

  REAL    :: deltaL(nLin,nLin),Lam(nLin,nLin),Lap(nLin,nLin) 

  REAL    :: RL(nVar,nLin),iRL(nLin,nVar),LLin(nLin,nLin) 

  REAL    :: RL(nVar,nLin),iRL(nLin,nVar), temp(nLin,nVar), LLin(nLin,nLin) 

  REAL    :: Aroe(nVar,nVar),Aroep(nVar,nVar), Aroem(nVar,nVar), Dm(nVar), Dp(nVar), dQ(nVar)

  REAL :: f1R(nVar), g1R(nVar), h1R(nVar) 

  REAL :: f1L(nVar), g1L(nVar), h1L(nVar) , VM(nVar) 

  REAL :: f1L(nVar), g1L(nVar), h1L(nVar) 

  !  

  nv(:)=0.

  nv(NormalNonZero+1)=1.

  !

  flattener=1.0 !0.8

  flattener=1.

  !

  CALL PDEFlux(f1L,g1L,h1L,QL)

  CALL PDEFlux(f1R,g1R,h1R,QR)
@@ -532,7 +502,7 @@ RECURSIVE SUBROUTINE HLLEMFluxFV(FL,FR,QL,QR,QavL,QavR,NormalNonZero) 
!ENDIF

  !USE Parameters, ONLY : d,nVar,ndim 

  QM = 0.5*(QL+QR) 

  CALL PDECons2PrimGRMHD(VM,QM,iErr)

  !CALL PDECons2PrimGRMHD(VM,QM,iErr)

  CALL PDEEigenvalues(LL,QL,nv)  

  CALL PDEEigenvalues(LR,QR,nv)  

!IF(ANY(QM(6:8).NE.0)) THEN
@@ -540,32 +510,24 @@ RECURSIVE SUBROUTINE HLLEMFluxFV(FL,FR,QL,QR,QavL,QavR,NormalNonZero) 
!    STOP

!ENDIF

  CALL PDEEigenvalues(LM,QM,nv)  

  sL = MIN( 0., MINVAL(LL(:)), MINVAL(LM(:)) ) 

  sR = MAX( 0., MAXVAL(LR(:)), MAXVAL(LM(:)) ) 

  minL = MINVAL(LL)

  minM = MINVAL(LM)

  maxR = MAXVAL(LR) 

  maxM = MAXVAL(LM) 

  sL = MIN( 0., minL, minM ) 

  sR = MAX( 0., maxR, maxM ) 

 ! PRINT *, "PDEIntermediateFields"

  !DO i=1,nVar

  !  WRITE(*,'(E16.6)'), QM(i)

  !ENDDO

  !print *,"*********************************************"

  !WRITE(*,'(a,f18.10,f18.10,f18.10)')    "***** nv:",nv(1),nv(2),nv(3)

  CALL PDEIntermediateFields(RL,LLin,iRL,QM,nv) 

  !PRINT *, "PDEIntermediateFields finished"

  Lam = 0.5*(LLin-ABS(LLin))

  Lap = 0.5*(LLin+ABS(LLin)) 

  deltaL = 0.0

  !print *,"*********************************************"

  !!print *,"*****LLin, QR(1),nv",QM(1),NormalNonZero

  !WRITE(*,'(a,E16.6,i9)')    "***** LLin, QR(1),nv",QM(1),NormalNonZero

  !print *,"**********************************************"



  DO i = 1, nLin

      deltaL(i,i) = (1. - Lam(i,i)/(sL-1e-14) - Lap(i,i)/(sR+1e-14) )*flattener(i)  

      !print *,"i,DeltaL(i,i):",i,DeltaL(i,i)

      !WRITE(*,'(a,i9,E16.6,E16.6,E16.6,E16.6,E16.6)') "i,QM(i),VM(i),DeltaL(i,i),LM(i):",i,QM(i),VM(i),DeltaL(i,i),LLin(i,i),LM(i)

      !WRITE(*,'(a, i9, f18.10, f16.5, f16.5, i9)') 

  ENDDO    

  !print *,"**********************************************"

  !STOP

#ifdef VISCOUS

  CALL PDEViscEigenvalues(LL,QL,nv)  

  CALL PDEViscEigenvalues(LR,QR,nv)
@@ -581,13 +543,10 @@ RECURSIVE SUBROUTINE HLLEMFluxFV(FL,FR,QL,QR,QavL,QavR,NormalNonZero) 
  IF(QR(1).LT.1e-9.OR.QL(1).LT.1e-9) THEN

      deltaL = 0.

  ENDIF

  temp = MATMUL( deltaL, iRL ) 

  absA = absA - sR*sL/(sR-sL)*MATMUL( RL, temp )  ! HLLEM anti-diffusion  

  !    

  !deltaL = 0.

  absA = absA - sR*sL/(sR-sL)*MATMUL( RL, MATMUL(deltaL, iRL) )  ! HLLEM anti-diffusion  

  !    

  !PRINT *, "RoeMatrix"

  CALL RoeMatrix(ARoe,QL,QR,nv)

  !PRINT *, "RoeMatrix done!"

  !

  ARoem = -sL*ARoe/(sR-sL)

  ARoep = +sR*ARoe/(sR-sL)

ApplicationExamples/GRMHDb/GRMHDb_ADERDG_FV_2D_clean/TecplotWriter.cpp

+0 −1
Original line number Diff line number Diff line
@@ -41,7 +41,6 @@ void GRMHDb::TecplotWriter::plotADERDGPatch( 
	plotForADERSolver = 1;

	plotForFVSolver = 0;

	elementcalltecplotaderdgplotter_(u, &offsetOfPatch[0], &sizeOfPatch[0], &plotForADERSolver);

	//printf("SONO QUI IN plotADERDGPatch");

}

ApplicationExamples/GRMHDb/GRMHDb_ADERDG_FV_2D_clean/TECPLOTPLOTTER.f90

+5 −5

File changed.

Contains only whitespace changes.

LRZ Homepage | Datenschutz | Dokumentation und Betriebsbedingungen | Impressum