Transformation of nodal value in RK44

Solver/CMakeLists.txt

@@ -23,6 +23,7 @@ set(SRC
   dgConservationLawMaxwell.cpp
   dgConservationLawPerfectGas.cpp
   dgLimiter.cpp
+  dgTransformNodalValue.cpp
   dgResidual.cpp
   dgDofContainer.cpp
   function.cpp

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)) 

Solver/dgAlgorithm.cpp

@@ -9,6 +9,7 @@
 #include "MEdge.h"
 #include "function.h"
 #include "dgLimiter.h"
+#include "dgTransformNodalValue.h"
 #include "meshPartition.h"
 
 /*

Solver/dgAlgorithm.h

@@ -11,6 +11,7 @@ class dgConservationLaw;
 class dgDofContainer;
 class dgTerm;
 class dgLimiter;
+class dgTransformNodalValue;
 class dgSystemOfEquations;
 
 class dgAlgorithm {

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");
+
 }
 

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);

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);

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();