diff --git a/Solver/TESTCASES/CylinderEddies.lua b/Solver/TESTCASES/CylinderEddies.lua
index cd73f5cf70b3b8a22126b1ed28576d79de2aa293..4969a88487f21408d4ff03ece0c7d52679c389c9 100644
--- a/Solver/TESTCASES/CylinderEddies.lua
+++ b/Solver/TESTCASES/CylinderEddies.lua
@@ -26,9 +26,9 @@ FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
-- diffusivity
mu=fullMatrix(1,1);
-mu:set(0,0,0.03)
+mu:set(0,0,0.02)
kappa=fullMatrix(1,1);
-kappa:set(0,0,0.03)
+kappa:set(0,0,0.02)
print'*** Loading the mesh and the model ***'
myModel = GModel ()
@@ -58,19 +58,20 @@ DG:exportSolution('output/cyl_0')
print'*** solve ***'
-CFL = 2.1;
+CFL = 20.1;
dt = CFL * DG:computeInvSpectralRadius();
print('DT = ',dt)
T = 0;
-for i=1,10000 do
+for i=1,100000 do
dt = CFL * DG:computeInvSpectralRadius();
norm = DG:RK44(dt)
T = T + dt
- if (i % 1 == 0) then
+ if (i % 10 == 0) then
print('*** ITER ***',i,norm,dt,T)
end
if (i % 100 == 0) then
DG:exportSolution(string.format("output/cyl-%06d", i))
+ DG:saveSolution(string.format("output/cyl-%06d.bin", i))
end
end
diff --git a/Solver/TESTCASES/ForwardFacingStep.lua b/Solver/TESTCASES/ForwardFacingStep.lua
index 196e2bab8ea92dec3b1dc9489928590a05b67b9a..d0ec53d98bd660d29ed094800943b47d89c5fb12 100644
--- a/Solver/TESTCASES/ForwardFacingStep.lua
+++ b/Solver/TESTCASES/ForwardFacingStep.lua
@@ -1,10 +1,11 @@
-MACH = 1.0;
+MACH = 3;
GAMMA = 1.4;
U = 3.0
V = 0.0
RHO = 1.4;
PRES = RHO*U*U/(GAMMA*MACH*MACH)
+
--PRES = 1;
--PRES = ./(MACH*RHO*RHO*GAMMA*GAMMA)
@@ -18,7 +19,7 @@ function free_stream( XYZ, FCT )
FCT:set(i,0,RHO)
FCT:set(i,1,RHO*U)
FCT:set(i,2,RHO*V)
- FCT:set(i,3, 0.5*RHO*(U*U+V*V)+PRES/(RHO*GAMMA-1))
+ FCT:set(i,3, 0.5*RHO*(U*U+V*V)+PRES/(GAMMA-1))
end
end
@@ -41,6 +42,7 @@ DG:setConservationLaw(law)
law:addBoundaryCondition('Walls',law:newWallBoundary())
law:addBoundaryCondition('LeftRight',law:newOutsideValueBoundary(FS))
+--law:addBoundaryCondition('Walls',law:newOutsideValueBoundary(FS))
DG:setup()
@@ -54,13 +56,17 @@ print'*** export ***'
DG:exportSolution('output/solution_0')
print'*** solve ***'
-CFL = 2.0;
-for i=1,1000 do
- dt = CFL * DG:computeInvSpectralRadius();
- norm = DG:RK44_limiter(0.1*dt)
- print('*** ITER ***',i,dt, norm)
- if (i % 1 == 0) then
+CFL = 2
+
+for i=1,5000 do
+ dt = CFL * DG:computeInvSpectralRadius();
+ norm = DG:RK44_limiter(dt)
+ DG:limitSolution()
+ if (i % 10 == 0) then
+ print('*** ITER ***',i,norm,dt)
+ end
+ if (i % 100 == 0) then
DG:exportSolution(string.format("output/solution-%06d", i))
end
end
diff --git a/Solver/TESTCASES/RayleighTaylor.lua b/Solver/TESTCASES/RayleighTaylor.lua
index f58f88d7762b500864b89a6165e864a1292a835b..dfd59a504b84722277ac6583f10d0aed9eb20757 100644
--- a/Solver/TESTCASES/RayleighTaylor.lua
+++ b/Solver/TESTCASES/RayleighTaylor.lua
@@ -39,7 +39,7 @@ end
Example of a lua program driving the DG code
--]]
-order = 4
+order = 1
print'*** Loading the mesh and the model ***'
myModel = GModel ()
myModel:load ('rect.geo')
diff --git a/Solver/TESTCASES/step.geo b/Solver/TESTCASES/step.geo
index 560bc22e2b48cafec8a9bdb0e6038f2dc4f095be..f7a3d05a96094d315e83a00aace6ed50d09510f1 100644
--- a/Solver/TESTCASES/step.geo
+++ b/Solver/TESTCASES/step.geo
@@ -1,8 +1,8 @@
Point(1) = {0, 0, 0, .1};
Point(2) = {.6, 0, 0, .1};
Point(3) = {.6, .2, 0, .1};
-Point(4) = {2.4, .2, 0, .1};
-Point(5) = {2.4, 1, 0, .1};
+Point(4) = {3, .2, 0, .1};
+Point(5) = {3, 1, 0, .1};
Point(6) = {0, 1, 0, .1};
Line(1) = {6, 5};
Line(2) = {5, 4};
@@ -13,5 +13,7 @@ Line(6) = {1, 6};
Line Loop(7) = {2, 3, 4, 5, 6, 1};
Plane Surface(8) = {7};
Physical Line("Walls") = {1, 3, 4, 5};
+//Physical Line("Walls") = {1};
Physical Line("LeftRight") = {2, 6};
+//Physical Line("LeftRight") = {2,3,4,5, 6};
Physical Surface("Body") = {8};
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 24a9de318ef0394b1b952c03ca66c3fd8403f57d..8d3dc1bac928221d85011ed5356f4c301e989e8c 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -281,6 +281,10 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
double b[4] = {0.,h/2.0,h/2.0,h};
+ if (orderRK = 1.0){
+ a[0] = h;
+ }
+
// fullMatrix<double> K(sol);
// Current updated solution
// fullMatrix<double> Unp(sol);
@@ -298,14 +302,17 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
Unp._data->axpy(*(sol._data));
for(int j=0; j<orderRK;j++){
- if(j){
- K._data->scale(b[j]);
- K._data->axpy(*(sol._data));
- if (limiter){
- limiter->apply(K, eGroups, fGroups);
+ if(j){
+ K._data->scale(b[j]);
+ K._data->axpy(*(sol._data));
+ if (limiter){
+ // limiter->apply(K, eGroups, fGroups);
}
- }
-
+ }
+
+ if (limiter){
+ limiter->apply(K, eGroups, fGroups);
+ }
this->residual(claw,eGroups,fGroups,bGroups,K._dataProxys,resd._dataProxys);
K._data->scale(0.);
for(int k=0;k < eGroups.size();k++) {
@@ -318,10 +325,10 @@ void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation l
}
}
Unp._data->axpy(*(K._data),a[j]);
+ //if (limiter) limiter->apply(Unp, eGroups, fGroups);
}
if (limiter) limiter->apply(Unp, eGroups, fGroups);
- for (int i=0;i<sol._dataSize;i++){
-
+ for (int i=0;i<sol._dataSize;i++){
(*sol._data)(i)=(*Unp._data)(i);
}
}
@@ -524,6 +531,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
normalFluxQP(iPt,k) = (*boundaryTerm)(iPt,k)*detJ;
}
+ // ----- 2.3.3 --- compute viscous contribution
if (diffusiveFluxLeft) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index b93a2777846e6b15074650471c6f7b1295816e6d..196b9b8b519a42a221573dc920edaa803163f041 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -43,6 +43,7 @@ class dgConservationLaw {
virtual dataCacheDouble *newConvectiveFlux (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newRiemannSolver (dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {return NULL;}
+ virtual dataCacheDouble *newClipToPhysics (dataCacheMap &cacheMap) const {return NULL;}
inline const dgBoundaryCondition *getBoundaryCondition(const std::string tag) const {
std::map<const std::string,dgBoundaryCondition*>::const_iterator it = _boundaryConditions.find(tag);
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 9160134f075f33120a58ad0fd4ee22540123b53a..d3cecab9427e3cf63998176841d4fa1ff8132e28 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -6,6 +6,206 @@
*/
const double GAMMA = 1.4;
+#define RIEM_TOL ((double)1e-12) /* Rel error to achieve in pstar */
+#define RIEM_ITER (20) /* Max number of iterations */
+
+#define smallu ((double)1e-14)
+#define smallp ((double)1e-14)
+
+#define max2(a,b) ( (a)>(b) ? (a) : (b) )
+#define max3(a,b,c) ( max2( a , max2(b,c) ) )
+
+bool riemannExact( double rho1, /* inputs */
+ double u1,
+ double v1,
+ double w1,
+ double p1,
+ double rho2,
+ double u2,
+ double v2,
+ double w2,
+ double p2,
+ double *rhoav, /* outputs */
+ double *uav,
+ double *utav,
+ double *uttav,
+ double *pav )
+{
+
+ /*
+ * local variables (previously in commons)
+ */
+ // double rho,u,p,ut,utt;
+
+ int n;
+ double plft, vlft, ulft, utlft, uttlft, clft;
+ double prght, vrght, urght, utrght, uttrght, crght;
+ double pstar, ustar, rhostr, cestar, westar, vstar;
+ double rhos, us, uts, utts, ps;
+ double ws, vs, wes, ces;
+ double wlft, wrght;
+ double scrch1, scrch2, scrch3, scrch4;
+ double prevpstar, relpstar;
+ // double Q[4];
+
+
+ //
+ // rml - convert to flash unknowns;
+ //
+
+ vlft= (double)1/ rho1; // volume over mass
+ ulft= u1;
+ utlft= v1;
+ uttlft= w1;
+ plft= p1;
+
+ if(plft/vlft < 0.0){
+ return false;}
+ clft=sqrt(GAMMA*plft/vlft); // rho * eulerian sound speed
+
+ vrght= (double)1/ rho2;
+ urght= u2;
+ utrght= v2;
+ uttrght= w2;
+ prght= p2;
+ if(prght/vrght < 0.0) {
+
+
+ return false;}
+ crght=sqrt(GAMMA*prght/vrght);
+
+ //
+ // Start of Bruce's original code
+ //
+
+
+ pstar=prght-plft-crght*(urght-ulft);
+ pstar=plft+pstar*clft/(clft+crght);
+ pstar=max2(smallp,pstar);
+
+ //
+ // rml - add ability to jump of out loop when converged
+ //
+
+ prevpstar=pstar; // save previous value
+ relpstar=RIEM_TOL;
+
+ for (n=0; n<RIEM_ITER && relpstar>=RIEM_TOL ; n++)
+ {
+ scrch1=0.5*(GAMMA+1.)/GAMMA;
+ wlft=1.+scrch1*(pstar-plft)/plft;
+ wrght=1.+scrch1*(pstar-prght)/prght;
+
+ if(wlft<0.0 || wrght<0.0)
+ {
+ return false;}
+ wlft=clft*sqrt(wlft);
+ wrght=crght*sqrt(wrght);
+
+ scrch1=4.*vlft*wlft*wlft;
+ scrch1=-scrch1*wlft/(scrch1-(GAMMA+1.)*(pstar-plft));
+ scrch2=4.*vrght*wrght*wrght;
+ scrch2=scrch2*wrght/(scrch2-(GAMMA+1.)*(pstar-prght));
+ scrch3=ulft-(pstar-plft)/wlft;
+ scrch4=urght+(pstar-prght)/wrght;
+ pstar=pstar+(scrch4-scrch3)*(scrch1*scrch2)/(scrch2-scrch1);
+ pstar=max2(smallp,pstar);
+
+ //
+ // rml - compute relative error
+ //
+
+ relpstar=fabs(pstar-prevpstar)/prevpstar;
+ // printf("pstar %e diff %e relpstar %e\n",
+ // pstar,pstar-prevpstar,relpstar);
+
+ prevpstar=pstar;
+ }
+
+ if (relpstar>=RIEM_TOL)
+ {
+ printf("Riemann solver failed : residual %12.5E\n",relpstar);
+ return false;
+ }
+
+ scrch3=ulft-(pstar-plft)/wlft;
+ scrch4=urght+(pstar-prght)/wrght;
+ ustar=0.5*(scrch3+scrch4);
+
+ if (ustar>=0)
+ scrch1= 1;
+ else
+ scrch1= -1;
+
+
+ scrch2=0.5*(1.+scrch1);
+ scrch3=0.5*(1.-scrch1);
+ ps=plft*scrch2+prght*scrch3;
+ us=ulft*scrch2+urght*scrch3;
+ uts=utlft*scrch2+utrght*scrch3;
+ utts=uttlft*scrch2+uttrght*scrch3; // rml 3rd dimension
+
+ vs=vlft*scrch2+vrght*scrch3;
+ rhos=1./vs;
+ ws=wlft*scrch2+wrght*scrch3;
+
+ if(ps*vs < 0) {
+ return false;
+ }
+
+ ces=sqrt(GAMMA*ps*vs);
+ vstar=vs-(pstar-ps)/(ws*ws);
+ vstar=max2(vstar,(GAMMA-1.)/(GAMMA+1.)*vs);
+ rhostr=1./vstar;
+
+ if(pstar*vstar < 0)
+ {
+ return false;}
+
+ cestar=sqrt(GAMMA*pstar*vstar);
+ wes=ces-scrch1*us;
+ westar=cestar-scrch1*ustar;
+ scrch4=ws*vs-scrch1*us;
+
+ if ( (pstar-ps) >= 0 )
+ {
+ wes=scrch4;
+ westar=scrch4;
+ }
+
+
+ //
+ // compute correct state for rarefaction fan
+ //
+
+ scrch1=max3(wes-westar,wes+westar,smallu);
+ scrch1=(wes+westar)/scrch1;
+ scrch1=0.5*(1.+scrch1);
+ scrch2=1.-scrch1;
+
+ *rhoav= scrch1*rhostr+scrch2*rhos;
+ *uav= scrch1*ustar+scrch2*us;
+ *utav= uts;
+ *uttav= utts; // rml 3rd dimension //
+ *pav=scrch1*pstar+scrch2*ps;
+
+ if (westar>=0)
+ {
+ *rhoav= rhostr;
+ *uav= ustar;
+ *pav= pstar;
+ }
+
+ if (wes<0)
+ {
+ *rhoav= rhos;
+ *uav= us;
+ *pav= ps;
+ }
+
+ return true;
+}
+
static inline void _ROE2D (const double &_GAMMA,
const double &nx,
const double &ny,
@@ -289,6 +489,72 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
}
};
+class dgPerfectGasLaw2d::riemannGodunov : public dataCacheDouble {
+ dataCacheDouble &normals, &solL, &solR;
+ public:
+ riemannGodunov(dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight):
+ normals(cacheMapLeft.get("Normals", this)),
+ solL(cacheMapLeft.get("Solution", this)),
+ solR(cacheMapRight.get("Solution", this))
+ {};
+ void _eval () {
+ int nQP = solL().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,12);
+
+ for(int i=0; i< nQP; i++) {
+ const double nx = normals(0,i);
+ const double ny = normals(1,i);
+
+ const double unL = (solL(i,1)*nx+solL(i,2)*ny)/solL(i,0);
+ const double utL = (solL(i,1)*ny-solL(i,2)*nx)/solL(i,0);
+ const double rhoV2L = (unL*unL+utL*utL)*solL(i,0);
+ const double pL = (GAMMA-1.)*solL(i,3) - 0.5*(GAMMA-1.)* rhoV2L;
+
+ const double unR = (solR(i,1)*nx+solR(i,2)*ny)/solR(i,0);
+ const double utR = (solR(i,1)*ny-solR(i,2)*nx)/solR(i,0);
+ const double rhoV2R = (unR*unR+utR*utR)*solR(i,0);
+ const double pR = (GAMMA-1.)*solR(i,3) - 0.5*(GAMMA-1.)* rhoV2R;
+
+ double rhoAv,unAv,utAv,usAv,pAv;
+ riemannExact(solL(i,0), unL,utL,0.0,pL,
+ solR(i,0), unR,utR,0.0,pR,
+ &rhoAv,&unAv,&utAv,&usAv,&pAv);
+
+
+ const double vxAv = unAv * nx + utAv *ny;
+ const double vyAv = unAv * ny - utAv *nx;
+ // p = rho E (G-1) - 0.5 (G-1) * rho V^2
+ // rho E = p / (G-1) + 0.5 * rho V^2
+ const double rhoE =
+ (1./(GAMMA-1.)) * pAv + 0.5 * rhoAv * (vxAv*vxAv + vyAv*vyAv);
+
+ /* printf("%g %g %g %g, (%g %g) %g %g %g %g\n",
+ solL(i,0), solL(i,1), solL(i,2), solL(i,3),pL,pAv,
+ rhoAv, rhoAv*vxAv,rhoAv*vyAv,rhoE);
+ */
+ const double F0 = rhoAv * unAv;
+ /*
+ const double qq = invrho*(sol(k,3)+p);
+
+ _value(k,0) = sol(k,1);
+ _value(k,1) = q11+p;
+ _value(k,2) = q12;
+ _value(k,3) = sol(k,1)*qq;
+ */
+
+ _value(i,0) = -F0;
+ _value(i,1) = -(F0*vxAv + pAv*nx);
+ _value(i,2) = -(F0*vyAv + pAv*ny);
+ _value(i,3) = - (rhoE+pAv) * unAv;
+ _value(i,4) = -_value(i,0);
+ _value(i,5) = -_value(i,1);
+ _value(i,6) = -_value(i,2);
+ _value(i,7) = -_value(i,3);
+ }
+ }
+};
+
class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
dataCacheDouble /
public:
@@ -341,7 +607,7 @@ dataCacheDouble *dgPerfectGasLaw2d::newConvectiveFlux( dataCacheMap &cacheMap) c
return new advection(cacheMap);
}
dataCacheDouble *dgPerfectGasLaw2d::newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
- return new riemann(cacheMapLeft, cacheMapRight);
+ return new riemannGodunov(cacheMapLeft, cacheMapRight);
}
dataCacheDouble *dgPerfectGasLaw2d::newDiffusiveFlux( dataCacheMap &cacheMap) const {
if (_muFunctionName.empty() || _kappaFunctionName.empty())
@@ -384,6 +650,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
for(int i=0; i< nQP; i++) {
const double nx = normals(0,i);
const double ny = normals(1,i);
+
const double solLeft [4] = {sol(i,0),sol(i,1),sol(i,2),sol(i,3)};
const double vn = (solLeft [1] * nx + solLeft [2] * ny);
const double solRight[4] = {sol(i,0),
@@ -399,12 +666,13 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
/*
const double q11 = sol(i,1)*sol(i,1)/sol(i,0);
const double q22 = sol(i,2)*sol(i,2)/sol(i,0);
- const double p = (GAMMA-1)*sol(i,3) - 0.5*(GAMMA-1)*(q11+q22);
+ const double p = (GAMMA-1.)*sol(i,3) - 0.5*(GAMMA-1.)*(q11+q22);
_value(i,0) = 0;//FLUX[0];
_value(i,1) = -p*nx;//FLUX[1];
_value(i,2) = -p*ny;//FLUX[2];
_value(i,3) = 0.0;//FLUX[3];
*/
+
}
}
};
@@ -474,7 +742,7 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
//-------------------------------------------------------------------------------
// Slip Wall Dirichlet BC --
-// Assume zero derivatives of all variables --> put the same as inside
+// Assume zero normal derivatives of all variables --> put the same as inside
//-------------------------------------------------------------------------------
class dirichletSlip : public dataCacheDouble {
diff --git a/Solver/dgConservationLawPerfectGas.h b/Solver/dgConservationLawPerfectGas.h
index 75e43d9a2c716868b1017da666c118780b462cbc..8eb66fd4b01e880027c209afbed128cf6696ad20 100644
--- a/Solver/dgConservationLawPerfectGas.h
+++ b/Solver/dgConservationLawPerfectGas.h
@@ -5,6 +5,7 @@ class dgPerfectGasLaw2d : public dgConservationLaw {
class advection;
class diffusion;
class riemann;
+ class riemannGodunov;
class source;
class maxConvectiveSpeed;
class maxDiffusivity;
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index 5225862323f3995216a6c3a1d01a85e525c8e97c..ba4f410d282cafc80c2a8e35ef1298612a0433d8 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -11,6 +11,7 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
//WARNING: ONLY FOR 1 GROUP OF FACES
//TODO: make this more general
+
dgGroupOfFaces* group = fGroups[0];
fullMatrix<double> &SolLeft = *(solution._dataProxys[0]);
fullMatrix<double> &SolRight = *(solution._dataProxys[0]);
@@ -18,19 +19,7 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
int totNbElems = solution.getNbElements();
- //--- CLIPPING: check unphysical values
- // double fmin = 1.e-10;
-// for (int iElement=0;iElement<totNbElems;iElement++){
-// fullMatrix<double> solElem;
-// solElem.setAsProxy(SolRight, nbFields*iElement, nbFields );
-// for (int k=0;k< solElem.size1() ;k++){
-// if (solElem(k,0) < fmin) {
-// printf("ARGG negative density \n");
-// solElem(k,0) = fmin;
-// }
-// }
-// }
-
+ // --- CLIPPING: check unphysical values
// first compute max and min of all fields for all stencils
//----------------------------------------------------------
@@ -39,65 +28,36 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
MIN.setAll ( 1.e22);
MAX.setAll (-1.e22);
- fullMatrix<double> MEANL (totNbElems, nbFields);
- fullMatrix<double> MEANR (totNbElems, nbFields);
- MEANL.setAll ( 0.01);
- MEANR.setAll ( 0.01);
-
int iElementL, iElementR, fSize;
for(int iFace=0 ; iFace<group->getNbElements();iFace++)
{
iElementL = group->getElementLeftId(iFace);
iElementR = group->getElementRightId(iFace);
- if (iElementR >= 0)
- {
- fullMatrix<double> TempL, TempR;
- TempL.setAsProxy(SolLeft, nbFields*iElementL, nbFields );
- TempR.setAsProxy(SolRight, nbFields*iElementR, nbFields );
-
- fSize = TempL.size1();
- for (int k=0; k< nbFields; ++k)
- {
- double AVGL = 0;
- double AVGR = 0;
- for (int i=0; i<fSize; ++i)
- {
- AVGL += TempL(i,k);
- AVGR += TempR(i,k);
- }
- AVGL /= (double) fSize;
- AVGR /= (double) fSize;
- MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
- MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
- MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
- MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
-
- MEANR( iElementL,k ) = AVGR;
- MEANL( iElementR,k ) = AVGL;
-
- }
- }
- else{
- fullMatrix<double> TempL ;
- TempL.setAsProxy(SolLeft, nbFields*iElementL, nbFields );
- for (int k=0; k<nbFields; ++k){
- for (int i=0; i<fSize; ++i){
- MIN ( iElementL , k ) = std::min ( TempL(i,k) , MIN ( iElementL , k ) );
- MAX ( iElementL , k ) = std::max ( TempL(i,k) , MAX ( iElementL , k ) );
-
-
- double AVG = 0;
- for (int i=0; i<fSize; ++i) AVG += TempL(i,k);
-
- }
- }
- }
- }
- // some parallel should be done here in order to compute averages on other processors
- //----------------------------------------------------------
- //...
+ fullMatrix<double> TempL, TempR;
+ TempL.setAsProxy(SolLeft, nbFields*iElementL, nbFields );
+ TempR.setAsProxy(SolRight, nbFields*iElementR, nbFields );
+
+ fSize = TempL.size1();
+ for (int k=0; k< nbFields; ++k){
+ double AVGL = 0;
+ double AVGR = 0;
+ for (int i=0; i<fSize; ++i) {
+ AVGL += TempL(i,k);
+ AVGR += TempR(i,k);
+ }
+ AVGL /= (double) fSize;
+ AVGR /= (double) fSize;
+ MIN ( iElementL , k ) = std::min ( AVGR , MIN ( iElementL , k ) );
+ MAX ( iElementL , k ) = std::max ( AVGR , MAX ( iElementL , k ) );
+ MIN ( iElementR , k ) = std::min ( AVGL , MIN ( iElementR , k ) );
+ MAX ( iElementR , k ) = std::max ( AVGL , MAX ( iElementR , k ) );
+ }
+ }
+ // printf("fSize = %d\n",fSize);
+
+ //----------------------------------------------------------
// then limit the solution
//----------------------------------------------------------
@@ -129,13 +89,42 @@ bool dgSlopeLimiter::apply ( dgDofContainer &solution,
if (locMin != AVG && locMin < neighMin) slopeLimiterValue = std::min ( slopeLimiterValue , (AVG-neighMin) / (AVG-locMin) );
if (AVG < neighMin) slopeLimiterValue = 0;
if (AVG > neighMax) slopeLimiterValue = 0;
+
+ // if (slopeLimiterValue != 1.0){
+ // printf("LIMTING %g\n",slopeLimiterValue);
+ // }
+ // slopeLimiterValue = 0.0;
- for (int i=0; i<fSize; ++i) Temp(i,k) *= slopeLimiterValue;
-
- for (int i=0; i<fSize; ++i) Temp(i,k) += AVG;
+ for (int i=0; i<fSize; ++i) Temp(i,k) = AVG + Temp(i,k)*slopeLimiterValue;
}
}
+
+ #if 0
+ double rhomin = 1.e-3;
+ double presmin= 1.e-3;
+ for (int iElement=0;iElement<totNbElems;iElement++){
+ fullMatrix<double> solElem;
+ solElem.setAsProxy(SolRight, nbFields*iElement, nbFields );
+ for (int k=0;k< solElem.size1() ;k++){
+ if (solElem(k,0) < rhomin) {
+ solElem(k,0) = rhomin;
+ }
+ double rhoV2 = 0;
+ for (int j=0;j<2;j++) {
+ double rhov = solElem(k,j+1);
+ rhoV2 += rhov*rhov;
+ }
+ rhoV2 /= solElem(k,0);
+ const double p = (1.4-1)*solElem(k,3) - 0.5*(1.4-1)* rhoV2;
+ if (p < presmin) {
+ solElem(k,3) = 0.5 *rhoV2 + presmin / (1.4-1);
+ // printf("negative pressure %g cliiped !!\n",p);
+ }
+ }
+ }
+ #endif
+
return true;
}
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index f8350962bb529f5a9f6d8eb34f32ce62ba74f36f..2d30b2923b6ae9dccf078263d6e787e47a048007 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -60,6 +60,7 @@ void dgSystemOfEquations::registerBindings(binding *b){
cm->setArgNames("functionName",NULL);
cm->setDescription("project the function \"functionName\" on the solution vector");
cm = cb->addMethod("RK44",&dgSystemOfEquations::RK44);
+ cm = cb->addMethod("ForwardEuler",&dgSystemOfEquations::ForwardEuler);
cm->setArgNames("norm","dt",NULL);
cm->setDescription("do a runge-kuta temporal iteration with a time step \"dt\" and return the sum of the nodal residuals");
cm = cb->addMethod("setOrder",&dgSystemOfEquations::setOrder);
@@ -114,11 +115,16 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
}
double dgSystemOfEquations::RK44_limiter(double dt){
- dgLimiter *sl = new dgSlopeLimiter();
- _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
- return _solution->_data->norm();
+ dgLimiter *sl = new dgSlopeLimiter();
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, sl);
+ delete sl;
+ return _solution->_data->norm();
}
+double dgSystemOfEquations::ForwardEuler(double dt){
+ _algo->rungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide, NULL,1);
+ return _solution->_data->norm();
+}
double dgSystemOfEquations::multirateRK43(double dt){
_algo->multirateRungeKutta(*_claw, _elementGroups, _faceGroups, _boundaryGroups, dt, *_solution, *_rightHandSide);
return _solution->_data->norm();
@@ -145,13 +151,13 @@ dgSystemOfEquations::~dgSystemOfEquations(){
}
}
-void dgSystemOfEquations::saveSolution (const std::string &name) const{
+void dgSystemOfEquations::saveSolution (std::string name) {
FILE *f = fopen (name.c_str(),"wb");
_solution->_data->binarySave(f);
fclose(f);
}
-void dgSystemOfEquations::loadSolution (const std::string &name){
+void dgSystemOfEquations::loadSolution (std::string name){
FILE *f = fopen (name.c_str(),"rb");
_solution->_data->binaryLoad(f);
fclose(f);
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 5799445d4d253de7cd3aa0ac7b15fa17e2a14282..22b11b795c8b68790a440308f8b9a70f9dc6bf78 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -60,6 +60,7 @@ public:
void limitSolution (); // apply the limiter on the solution
double RK44 (double dt);
double RK44_limiter (double dt);
+ double ForwardEuler (double dt);
double multirateRK43 (double dt);
void L2Projection (std::string functionName); // assign the solution to a given function
double computeInvSpectralRadius();
@@ -71,11 +72,10 @@ public:
iElement * _claw->nbFields(),_claw->nbFields());
}
void export_solution_as_is (const std::string &fileName);
- void saveSolution (const std::string &fileName) const;
- void loadSolution (const std::string &fileName);
+ void saveSolution (std::string fileName) ;
+ void loadSolution (std::string fileName);
~dgSystemOfEquations();
private:
};
-
#endif // _DG_SYSTEM_OF_EQUATIONS_