diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 6aa81f0fe1bcf54dce1049a6311ad3125b027bc4..f5c4fc47a2d79cb34aee0137676000e28a7056de 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -23,6 +23,7 @@ set(SRC
dgConservationLawMaxwell.cpp
dgConservationLawPerfectGas.cpp
dgLimiter.cpp
+ dgTransformNodalValue.cpp
dgResidual.cpp
dgDofContainer.cpp
function.cpp
diff --git a/Solver/TESTCASES/Diffusion.lua b/Solver/TESTCASES/Diffusion.lua
index b8093035934f6b9b44b08ab2bc77f13ffefe114b..2b2800dcf5df3189a6805c3a6ff1e32d35567f32 100644
--- a/Solver/TESTCASES/Diffusion.lua
+++ b/Solver/TESTCASES/Diffusion.lua
@@ -2,14 +2,14 @@ model = GModel ()
model:load ('square.geo')
model:load ('square.msh')
dg = dgSystemOfEquations (model)
-dg:setOrder(3)
+dg:setOrder(1)
-- conservation law
-- advection speed
nu=fullMatrix(1,1);
nu:set(0,0,0.01)
-law = dgConservationLawAdvection('',functionConstant(nu):getName())
+law = dgConservationLawAdvectionDiffusion('',functionConstant(nu):getName())
dg:setConservationLaw(law)
-- boundary condition
@@ -34,7 +34,8 @@ dg:exportSolution('output/Diffusion_00000')
-- main loop
for i=1,10000 do
- norm = dg:RK44(0.001)
+ --norm = dg:RK44(0.001)
+ norm = dg:RK44_TransformNodalValue(0.00001)
if (i % 50 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Diffusion-%05d", i))
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index a5a672144b95f369cf37419c7851e168f714f12f..8a833816711e3e67279d75227830c0d2a7c03551 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -9,6 +9,7 @@
#include "MEdge.h"
#include "function.h"
#include "dgLimiter.h"
+#include "dgTransformNodalValue.h"
#include "meshPartition.h"
/*
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 5756cce915bcae4ddfd70b99e682788307d1b9d7..e5bb4c943471f9096c3215755d4c264667c10b05 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -11,6 +11,7 @@ class dgConservationLaw;
class dgDofContainer;
class dgTerm;
class dgLimiter;
+class dgTransformNodalValue;
class dgSystemOfEquations;
class dgAlgorithm {
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
index 4f22873c08dac89f8dcf60f8232be15487268548..ffcf59692b78a40d116f8ec05d0262a056704622 100644
--- a/Solver/dgRungeKutta.cpp
+++ b/Solver/dgRungeKutta.cpp
@@ -3,6 +3,7 @@
#include "dgConservationLaw.h"
#include "dgDofContainer.h"
#include "dgLimiter.h"
+#include "dgTransformNodalvalue.h"
#include "dgResidual.h"
#include "dgGroupOfElements.h"
@@ -108,10 +109,14 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
dgDofContainer K (groups, nbFields);
dgDofContainer Unp (groups, nbFields);
dgDofContainer resd(groups, nbFields);
- K.scale(0.);
- K.axpy(*sol);
Unp.scale(0.);
Unp.axpy(*sol);
+
+ if (_TransformNodalValue) _TransformNodalValue->apply(&Unp);
+
+ K.scale(0.);
+ K.axpy(Unp);
+
dgResidual residual(claw);
for(int j=0; j<nStages;j++){
@@ -119,7 +124,7 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
K.scale(A[j]*dt);
K.axpy(*sol);
}
-
+ if (_TransformNodalValue) _TransformNodalValue->apply_Inverse(&K);
if (_limiter) _limiter->apply(&K);
residual.compute(*groups,K,resd);
K.scale(0.);
@@ -136,6 +141,7 @@ double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
Unp.axpy(K,b[j]*dt);
}
if (_limiter) _limiter->apply(&Unp);
+ if (_TransformNodalValue) _TransformNodalValue->apply_Inverse(&Unp);
sol->scale(0.);
sol->axpy(Unp);
return sol->norm();
@@ -241,5 +247,10 @@ void dgRungeKutta::registerBindings(binding *b) {
cm = cb->addMethod("setLimiter",&dgRungeKutta::setLimiter);
cm->setArgNames("limiter",NULL);
cm->setDescription("if a limiter is set, it is applied after each RK stage");
+
+ cm = cb->addMethod("setTransformNodalValue",&dgRungeKutta::setTransformNodalValue);
+ cm->setArgNames("TransformNodalValue",NULL);
+ cm->setDescription("if the Nodal values is transformed in first step of RK");
+
}
diff --git a/Solver/dgRungeKutta.h b/Solver/dgRungeKutta.h
index 6c9ec32a9cb93aeeb9f0c32858036eb0ace01899..8088e783c80f3916f85f3244302ac2a582bb1eb2 100644
--- a/Solver/dgRungeKutta.h
+++ b/Solver/dgRungeKutta.h
@@ -5,6 +5,7 @@
class dgConservationLaw;
class dgDofContainer;
class dgLimiter;
+class dgTransformNodalValue;
class dgGroupCollection;
class dgGroupOfElements;
class dgGroupOfFaces;
@@ -17,8 +18,11 @@ class dgRungeKutta {
double diagonalRK(const dgConservationLaw *claw, double dt, dgDofContainer *solution, int nStages, double *A, double *b); // c == A
double nonDiagonalRK(const dgConservationLaw *claw, double dt, dgDofContainer *solution, int nStages, fullMatrix<double>&A, double *b, double *c);
dgLimiter *_limiter;
+ dgTransformNodalValue *_TransformNodalValue;
public:
void setLimiter(dgLimiter *limiter) { _limiter = limiter; }
+ void setTransformNodalValue(dgTransformNodalValue *TransformNodalValue) { _TransformNodalValue = TransformNodalValue; }
+
double iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate22(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
double iterate33(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index a897fd46e039efe0cf31661b43ab014b1c412b64..7f7d62219c6a45791108ae0ffed1b65ecad3a78e 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -8,6 +8,7 @@
#include "PView.h"
#include "PViewData.h"
#include "dgLimiter.h"
+#include "dgTransformNodalValue.h"
#include "dgRungeKutta.h"
#ifdef HAVE_MPI
#include "mpi.h"
@@ -69,6 +70,16 @@ void dgSystemOfEquations::registerBindings(binding *b) {
cm = cb->addMethod("RK44_limiter",&dgSystemOfEquations::RK44_limiter);
cm->setArgNames("dt",NULL);
cm->setDescription("do one RK44 time step with the slope limiter (only for p=1)");
+
+cm = cb->addMethod("RK44_TransformNodalValue",&dgSystemOfEquations::RK44_TransformNodalValue);
+ cm->setArgNames("dt",NULL);
+ cm->setDescription("do one RK44 time step with a transformation of nodal value");
+
+ cm = cb->addMethod("RK44_TransformNodalValue_limiter",&dgSystemOfEquations::RK44_TransformNodalValue_limiter);
+ cm->setArgNames("dt",NULL);
+ cm->setDescription("do one RK44 time step with a limiter and a transformation of nodal value");
+
+
////cm = cb->addMethod("multirateRK43",&dgSystemOfEquations::multirateRK43);
//cm->setArgNames("dt",NULL);
//cm->setDescription("Do a runge-kuta temporal iteration with 4 sub-steps and a precision order of 3 using different time-step depending on the element size. This function returns the sum of the nodal residuals.");
@@ -126,6 +137,31 @@ double dgSystemOfEquations::RK44_limiter(double dt){
return _solution->norm();
}
+
+double dgSystemOfEquations::RK44_TransformNodalValue(double dt){
+ dgTransformNodalValue *tnv = new dgSupraTransformNodalValue(_claw);
+ dgRungeKutta rk;
+ rk.setTransformNodalValue(tnv);
+ rk.iterate44(_claw, dt, _solution);
+ delete tnv;
+ return _solution->norm();
+}
+
+
+
+double dgSystemOfEquations::RK44_TransformNodalValue_limiter(double dt){
+ dgLimiter *sl = new dgSlopeLimiter(_claw);
+ dgTransformNodalValue *tnv = new dgSupraTransformNodalValue(_claw);
+ dgRungeKutta rk;
+ rk.setTransformNodalValue(tnv);
+ rk.setLimiter(sl);
+ rk.iterate44(_claw, dt, _solution);
+ delete sl;
+ delete tnv;
+ return _solution->norm();
+}
+
+
double dgSystemOfEquations::ForwardEuler(double dt){
dgRungeKutta rk;
rk.iterateEuler(_claw, dt, _solution);
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 391402b9a5014cc7c1b64a696ef6b0b84a95468e..abc1a92abc79bff359b69b0923e00b5899d4ef0b 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -8,6 +8,7 @@
#include "dgConservationLaw.h"
#include "Gmsh.h"
#include "dgLimiter.h"
+#include "dgTransformNodalValue.h"
#include "dgDofContainer.h"
@@ -43,6 +44,9 @@ public:
void L2Projection (std::string functionName); // assign the solution to a given function
double RK44(double dt);
double RK44_limiter(double dt);
+ double RK44_TransformNodalValue(double dt);
+ double RK44_TransformNodalValue_limiter(double dt);
+
double ForwardEuler(double dt);
double computeInvSpectralRadius();