diff --git a/Solver/TESTCASES/ConservationTestMRAdvection.lua b/Solver/TESTCASES/ConservationTestMRAdvection.lua
index 734b937a0ef517c01ac50e40eb9a4fc0002ffa36..32cb24bd49afd7c1fa2a679689d672dd6e2d62e9 100644
--- a/Solver/TESTCASES/ConservationTestMRAdvection.lua
+++ b/Solver/TESTCASES/ConservationTestMRAdvection.lua
@@ -1,8 +1,6 @@
PI = 3.14159
+-- Defining functions
---[[
- Function for initial conditions
---]]
function initial_condition( XYZ, FCT )
for i=0,XYZ:size1()-1 do
X = XYZ:get(i,0)
@@ -11,6 +9,7 @@ function initial_condition( XYZ, FCT )
FCT:set(i,0,V)
end
end
+
function velocity( XYZ, FCT )
for i=0,XYZ:size1()-1 do
X = XYZ:get(i,0)
@@ -27,11 +26,9 @@ function toIntegrate( sol, FCT )
FCT:set(i,0,s)
end
end
---[[
- Example of a lua program driving the DG code
---]]
order = 1
+dimension = 2
print'*** Loading the mesh and the model ***'
myModel = GModel ()
myModel:load ('rect.geo')
@@ -48,23 +45,47 @@ elseif (order == 5) then
end
print'*** Create a dg solver ***'
+
+--functions
xyzF = functionCoordinates.get()
velF = functionLua(3, 'velocity', {xyzF})
nuF = functionConstant({0})
-law=dgConservationLawAdvectionDiffusion(velF,nuF)
+--Conservation laws
+law=dgConservationLawAdvectionDiffusion(velF,nuF)
+--Add the boundary conditions
law:addBoundaryCondition('Walls',law:new0FluxBoundary())
law:addBoundaryCondition('Top',law:new0FluxBoundary())
+
+--Fix initial condition
FS = functionLua(1, 'initial_condition', {xyzF})
+-- Creating the different "base" Element groups before splitting them (if multirate)
+-- Singlerate
+GC0=dgGroupCollection(myModel,dimension,order)
+solTmp0=dgDofContainer(GC0,1)
+solTmp0:L2Projection(FS)
+GC0:buildGroupsOfInterfaces();
+GC1=dgGroupCollection(myModel,dimension,order)
+solTmp1=dgDofContainer(GC1,1)
+solTmp1:L2Projection(FS)
+GC1:buildGroupsOfInterfaces();
-GC=dgGroupCollection(myModel,2,order)
-solTmp=dgDofContainer(GC,1)
-solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(2,1,2,law,solTmp)
+--Multirate
+GC2=dgGroupCollection(myModel,dimension,order)
+solTmp2=dgDofContainer(GC2,1)
+solTmp2:L2Projection(FS)
+
+GC3=dgGroupCollection(myModel,dimension,order)
+solTmp3=dgDofContainer(GC3,1)
+solTmp3:L2Projection(FS)
+
+GC4=dgGroupCollection(myModel,dimension,order)
+solTmp4=dgDofContainer(GC4,1)
+solTmp4:L2Projection(FS)
-- Example of Conservative Multirate scheme with Butcher tables as input (here RK2a)
A=fullMatrix(2, 2)
@@ -83,13 +104,14 @@ c:set(0, 1, 1)
-- Building the RK objects
-multirateRK=dgRungeKuttaMultirateConservative.new2b(GC, law)
-multirateRK:printButcher()
-- RK2a
+print'--- RK2a ---'
RK2a=dgRungeKutta()
-
+--dt0=RK2a:splitForMultirate(100, solTmp)
-- RK43
+print'--- RK43 ---'
RK43=dgRungeKutta()
+--dt1=RK43:splitForMultirate(100, solTmp)
-- Multirate RK2a
@@ -97,27 +119,37 @@ RK43=dgRungeKutta()
--multirateRK2a=dgRungeKuttaMultirate22(GC,law)
-- *** Multirate RK2a with as input the butcher tables ***
-multirateRK2a=dgRungeKuttaMultirateConservative(GC, law, A, b, c)
+print'--- multirateRK2a ---'
+multirateRK2a=dgRungeKuttaMultirateConservative(GC2, law, A, b, c)
+dt2=multirateRK2a:splitForMultirate(10, solTmp2)
-- Multirate RK43
-multirateRK43=dgRungeKuttaMultirate43(GC,law)
+print'--- multirateRK43 ---'
+multirateRK43=dgRungeKuttaMultirate43(GC3,law)
+dt3=multirateRK43:splitForMultirate(10, solTmp3)
-- Multirate RK2b
-multirateRK2b=dgRungeKuttaMultirateConservative.new2b(GC, law)
+print'--- multirateRK2b ---'
+multirateRK2b=dgRungeKuttaMultirateConservative.new2b(GC4, law)
-- multirateRK2b:printButcher()
+dt4=multirateRK2b:splitForMultirate(10, solTmp4)
+
+dt0=dt2
+dt1=dt2
+
+solution0=dgDofContainer(GC0,1)
+solution1=dgDofContainer(GC1,1)
+solution2=dgDofContainer(GC2,1)
+solution3=dgDofContainer(GC3,1)
+solution4=dgDofContainer(GC4,1)
-GC:buildGroupsOfInterfaces(myModel,2,order)
-solution0=dgDofContainer(GC,1)
-solution1=dgDofContainer(GC,1)
-solution2=dgDofContainer(GC,1)
-solution3=dgDofContainer(GC,1)
-solution4=dgDofContainer(GC,1)
solution0:L2Projection(FS)
solution1:L2Projection(FS)
solution2:L2Projection(FS)
solution3:L2Projection(FS)
solution4:L2Projection(FS)
-solution0:exportGroupIdMsh()
+
+solution4:exportGroupIdMsh()
print'---- Groups splitted ----'
@@ -131,8 +163,7 @@ solution4:exportMsh(string.format("outputMultirate/adv_mr2b-%06d", 0))
print'*** solve ***'
LC = 0.1*.1
-dt=dt
-print('DT=',dt)
+print('DT=',dt0)
norm0=solution0:norm()
norm1=solution1:norm()
@@ -141,42 +172,58 @@ norm3=solution3:norm()
norm4=solution4:norm()
print('Norm: ', norm0, norm1, norm2, norm3, norm4)
-dt=dt
-time=0
+time0=0
+time1=0
+time2=0
+time3=0
+time4=0
i=0
integrator=dgFunctionIntegrator(functionLua(1, 'toIntegrate', {functionSolution.get()}))
- print('*** ITER-START ***')
- w0=fullMatrix(1,1);
- w1=fullMatrix(1,1);
- w2=fullMatrix(1,1);
- w3=fullMatrix(1,1);
- w4=fullMatrix(1,1);
- integrator:compute(solution0,w0)
- integrator:compute(solution1,w1)
- integrator:compute(solution2,w2)
- integrator:compute(solution3,w3)
- integrator:compute(solution4,w4)
- w0ref=w0:get(0,0,0)
- w1ref=w1:get(0,0,0)
- w2ref=w2:get(0,0,0)
- w3ref=w3:get(0,0,0)
- w4ref=w4:get(0,0,0)
- print'*****************************************'
- print('integral (SR2a) : ',w0:get(0,0,0))
- print('integral (SR43) : ',w1:get(0,0,0))
- print('integral (MR2a) : ',w2:get(0,0,0))
- print('integral (MR43) : ',w3:get(0,0,0))
- print('integral (MR2b) : ',w4:get(0,0,0))
- print'******************************************'
-while i<10 do
- norm0 = RK2a:iterate22(law, dt, solution0)
- norm1 = RK43:iterate43SchlegelJCAM2009(law, dt, solution1)
- norm2 = multirateRK2a:iterate(dt, solution2)
- norm3 = multirateRK43:iterate(dt, solution3)
- norm4 = multirateRK2b:iterate(dt, solution4)
- time=time+dt
- if (i % 1 == 0) then
- print('*** ITER ***',i)--,time,norm0, norm1,norm2, norm3)
+
+print('*** ITER-START ***')
+w0=fullMatrix(1,1);
+w1=fullMatrix(1,1);
+w2=fullMatrix(1,1);
+w3=fullMatrix(1,1);
+w4=fullMatrix(1,1);
+
+integrator:compute(solution0,w0)
+integrator:compute(solution1,w1)
+integrator:compute(solution2,w2)
+integrator:compute(solution3,w3)
+integrator:compute(solution4,w4)
+
+w0ref=w0:get(0,0,0)
+w1ref=w1:get(0,0,0)
+w2ref=w2:get(0,0,0)
+w3ref=w3:get(0,0,0)
+w4ref=w4:get(0,0,0)
+
+print'*****************************************'
+print('integral (SR2a) : ',w0:get(0,0,0))
+print('integral (SR43) : ',w1:get(0,0,0))
+print('integral (MR2a) : ',w2:get(0,0,0))
+print('integral (MR43) : ',w3:get(0,0,0))
+print('integral (MR2b) : ',w4:get(0,0,0))
+print'******************************************'
+
+while i<100 do
+ --iterate
+ norm0 = RK2a:iterate22(law, dt0, solution0)
+ norm1 = RK43:iterate43SchlegelJCAM2009(law, dt1, solution1)
+ norm2 = multirateRK2a:iterate(dt2, solution2)
+ norm3 = multirateRK43:iterate(dt3, solution3)
+ norm4 = multirateRK2b:iterate(dt4, solution4)
+
+ --updating time
+ time0=time0+dt0
+ time1=time1+dt1
+ time2=time2+dt2
+ time3=time3+dt3
+ time4=time4+dt4
+
+ if (i % 10 == 0) then
+ print('*** ITER ***',i)
print'-----------------------------------------'
integrator:compute(solution0, w0)
integrator:compute(solution1, w1)
diff --git a/Solver/dgRungeKuttaMultirate.cpp b/Solver/dgRungeKuttaMultirate.cpp
index 65e5142fe378e0a461fc156eb7af510adf3a21c2..5c2392f39688ae3041d258e7723a5f6e4c870eec 100644
--- a/Solver/dgRungeKuttaMultirate.cpp
+++ b/Solver/dgRungeKuttaMultirate.cpp
@@ -285,6 +285,11 @@ void dgRungeKuttaMultirate::registerBindings(binding *b) {
dgRungeKuttaMultirate43::dgRungeKuttaMultirate43(dgGroupCollection *gc,dgConservationLaw* law): dgRungeKuttaMultirate(gc,law){}
+double dgRungeKuttaMultirate43::splitForMultirate(int maxLevels, dgDofContainer *solution){
+ double dt = _gc->splitGroupsForMultirate(maxLevels, 1, 1, _law, solution);
+ _gc->buildGroupsOfInterfaces();
+ return dt;
+}
double dgRungeKuttaMultirate43::iterate(double dt, dgDofContainer *solution){
initialize(10);
_solution=solution;
@@ -465,6 +470,9 @@ void dgRungeKuttaMultirate43::registerBindings(binding *b) {
cm = cb->setConstructor<dgRungeKuttaMultirate43,dgGroupCollection *,dgConservationLaw*>();
cm->setArgNames("groupCollection","law",NULL);
cm->setDescription("A new multirate explicit runge kutta, pass parameters to the iterate function");
+ cm = cb->addMethod("splitForMultirate",&dgRungeKuttaMultirate43::splitForMultirate);
+ cm->setArgNames("maxLevels","solution",NULL);
+ cm->setDescription("Split Groups for Multirate based on CFL condition, with appropriate buffer sizes, and build the corresponding Groups of Interfaces");
cm = cb->addMethod("iterate",&dgRungeKuttaMultirate43::iterate);
cm->setArgNames("dt","solution",NULL);
cm->setDescription("update the solution by doing a multirate RK43 (from Schlegel et al. Journal of Computational and Applied Mathematics, 2009) step of base time step dt for the conservation law");
@@ -473,6 +481,11 @@ void dgRungeKuttaMultirate43::registerBindings(binding *b) {
dgRungeKuttaMultirate22::dgRungeKuttaMultirate22(dgGroupCollection *gc,dgConservationLaw* law): dgRungeKuttaMultirate(gc,law){}
+double dgRungeKuttaMultirate22::splitForMultirate(int maxLevels, dgDofContainer *solution){
+ double dt = _gc->splitGroupsForMultirate(maxLevels, 1, 2, _law, solution);
+ _gc->buildGroupsOfInterfaces();
+ return dt;
+}
double dgRungeKuttaMultirate22::iterate(double dt, dgDofContainer *solution){
initialize(4);
_solution=solution;
@@ -627,6 +640,9 @@ void dgRungeKuttaMultirate22::registerBindings(binding *b) {
cm = cb->setConstructor<dgRungeKuttaMultirate22,dgGroupCollection *,dgConservationLaw*>();
cm->setArgNames("groupCollection","law",NULL);
cm->setDescription("A new multirate explicit runge kutta, pass parameters to the iterate function");
+ cm = cb->addMethod("splitForMultirate",&dgRungeKuttaMultirate22::splitForMultirate);
+ cm->setArgNames("maxLevels","solution",NULL);
+ cm->setDescription("Split Groups for Multirate based on CFL condition, with appropriate buffer sizes, and build the corresponding Groups of Interfaces");
cm = cb->addMethod("iterate",&dgRungeKuttaMultirate22::iterate);
cm->setArgNames("dt","solution",NULL);
cm->setDescription("update the solution by doing a multirate RK2a (from Constantinescu and Sandu, 'Update on Multirate Timestepping Methods for Hyperbolic Conservation Laws', Computer Sciance Technical Report, 2007) step of base time step dt for the conservation law");
@@ -686,6 +702,13 @@ void dgRungeKuttaMultirateConservative::printButcher(){
}
+
+double dgRungeKuttaMultirateConservative::splitForMultirate(int maxLevels, dgDofContainer *solution){
+ double dt = _gc->splitGroupsForMultirate(maxLevels, 1, _b->size2(), _law, solution);
+ _gc->buildGroupsOfInterfaces();
+ return dt;
+}
+
double dgRungeKuttaMultirateConservative::iterate(double dt, dgDofContainer *solution){
initialize(4);
_solution=solution;
@@ -937,22 +960,25 @@ void dgRungeKuttaMultirateConservative::registerBindings(binding *b) {
cm = cb->setConstructor<dgRungeKuttaMultirateConservative,dgGroupCollection*,dgConservationLaw* , fullMatrix<double>*, fullMatrix<double>*, fullMatrix<double>* >();
cm->setArgNames("groupCollection","law","A","b","c",NULL);
cm->setDescription("A new multirate explicit runge kutta, pass parameters to the iterate function");
+ cm = cb->addMethod("splitForMultirate",&dgRungeKuttaMultirateConservative::splitForMultirate);
+ cm->setArgNames("maxLevels","solution",NULL);
+ cm->setDescription("Split Groups for Multirate based on CFL condition, with appropriate buffer sizes, and build the corresponding Groups of Interfaces");
cm = cb->addMethod("iterate",&dgRungeKuttaMultirateConservative::iterate);
cm->setArgNames("dt","solution",NULL);
cm->setDescription("update the solution by doing a multirate RK constructed from the given Butcher tables (from Constantinescu and Sandu, 'Update on Multirate Timestepping Methods for Hyperbolic Conservation Laws', Computer Sciance Technical Report, 2007) step of base time step dt for the conservation law");
- cm = cb->addMethod("new44", &dgRungeKuttaMultirateConservative::new43);
+ cm = cb->addMethod("new44", &dgRungeKuttaMultirateConservative::new43);
cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK44");
cm->setArgNames("groupCollection", "law", NULL);
- cm = cb->addMethod("new43", &dgRungeKuttaMultirateConservative::new43);
+ cm = cb->addMethod("new43", &dgRungeKuttaMultirateConservative::new43);
cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK43");
cm->setArgNames("groupCollection", "law", NULL);
- cm = cb->addMethod("new2a", &dgRungeKuttaMultirateConservative::new2a);
+ cm = cb->addMethod("new2a", &dgRungeKuttaMultirateConservative::new2a);
cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK2a");
cm->setArgNames("groupCollection", "law", NULL);
- cm = cb->addMethod("new2b", &dgRungeKuttaMultirateConservative::new2b);
+ cm = cb->addMethod("new2b", &dgRungeKuttaMultirateConservative::new2b);
cm->setDescription("Creates a new Conservative Runge-Kutta scheme based on the base method RK2b");
cm->setArgNames("groupCollection", "law", NULL);
- cm = cb->addMethod("printButcher", &dgRungeKuttaMultirateConservative::printButcher);
+ cm = cb->addMethod("printButcher", &dgRungeKuttaMultirateConservative::printButcher);
cm->setDescription("Print the Butcher Tables for Buffers");
cm->setArgNames(NULL);
}
diff --git a/Solver/dgRungeKuttaMultirate.h b/Solver/dgRungeKuttaMultirate.h
index 100cac82032fa6f297339a37bbb5c5fcb70ea7ff..9aadec32caada764abbb4cb28f525283d2579111 100644
--- a/Solver/dgRungeKuttaMultirate.h
+++ b/Solver/dgRungeKuttaMultirate.h
@@ -25,6 +25,7 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
dgRungeKuttaMultirate(dgGroupCollection *gc,dgConservationLaw* law);
virtual ~dgRungeKuttaMultirate();
virtual double iterate(double dt, dgDofContainer *solution)=0;
+ virtual double splitForMultirate(int maxLevels, dgDofContainer *solution)=0;
static void registerBindings(binding *b);
};
@@ -32,6 +33,7 @@ class dgRungeKuttaMultirate43: public dgRungeKuttaMultirate{
void computeInputForK(int iK,int exponent,bool isBuffer);
void updateSolution(int exponent,bool isBuffer);
public:
+ double splitForMultirate(int maxLevels, dgDofContainer *solution);
double iterate(double dt, dgDofContainer *solution);
dgRungeKuttaMultirate43(dgGroupCollection *gc,dgConservationLaw *law);
static void registerBindings(binding *b);
@@ -41,31 +43,33 @@ class dgRungeKuttaMultirate22: public dgRungeKuttaMultirate{
void computeInputForK(int iK,int exponent,bool isBuffer,int upperLeveliK);
void updateSolution(int exponent,bool isBuffer);
public:
+ double splitForMultirate(int maxLevels, dgDofContainer *solution);
double iterate(double dt, dgDofContainer *solution);
dgRungeKuttaMultirate22(dgGroupCollection *gc,dgConservationLaw *law);
static void registerBindings(binding *b);
};
class dgRungeKuttaMultirateConservative: public dgRungeKuttaMultirate{
- private:
- fullMatrix<double> *_A;
- fullMatrix<double> *_b;
- fullMatrix<double> *_c;
- fullMatrix<double> *_AOuter;
- fullMatrix<double> *_bOuter;
- fullMatrix<double> *_cOuter;
- fullMatrix<double> *_AInner;
- fullMatrix<double> *_bInner;
- fullMatrix<double> *_cInner;
+ private:
+ fullMatrix<double> *_A;
+ fullMatrix<double> *_b;
+ fullMatrix<double> *_c;
+ fullMatrix<double> *_AOuter;
+ fullMatrix<double> *_bOuter;
+ fullMatrix<double> *_cOuter;
+ fullMatrix<double> *_AInner;
+ fullMatrix<double> *_bInner;
+ fullMatrix<double> *_cInner;
- void computeInputForK(int iK,int exponent,bool isBuffer,int upperLeveliK);
- void updateSolution(int exponent,bool isBuffer);
- public:
- double iterate(double dt, dgDofContainer *solution);
- dgRungeKuttaMultirateConservative(dgGroupCollection *gc,dgConservationLaw *law, fullMatrix<double> *A, fullMatrix<double> *b,fullMatrix<double> *c);
- void printButcher();
- static dgRungeKuttaMultirateConservative* new44(dgGroupCollection *gc, dgConservationLaw *law);
- static dgRungeKuttaMultirateConservative* new43(dgGroupCollection *gc, dgConservationLaw *law);
- static dgRungeKuttaMultirateConservative* new2a(dgGroupCollection *gc, dgConservationLaw *law);
- static dgRungeKuttaMultirateConservative* new2b(dgGroupCollection *gc, dgConservationLaw *law);
- static void registerBindings(binding *b);
+ void computeInputForK(int iK,int exponent,bool isBuffer,int upperLeveliK);
+ void updateSolution(int exponent,bool isBuffer);
+ public:
+ double splitForMultirate(int maxLevels, dgDofContainer *solution);
+ double iterate(double dt, dgDofContainer *solution);
+ dgRungeKuttaMultirateConservative(dgGroupCollection *gc,dgConservationLaw *law, fullMatrix<double> *A, fullMatrix<double> *b,fullMatrix<double> *c);
+ void printButcher();
+ static dgRungeKuttaMultirateConservative* new44(dgGroupCollection *gc, dgConservationLaw *law);
+ static dgRungeKuttaMultirateConservative* new43(dgGroupCollection *gc, dgConservationLaw *law);
+ static dgRungeKuttaMultirateConservative* new2a(dgGroupCollection *gc, dgConservationLaw *law);
+ static dgRungeKuttaMultirateConservative* new2b(dgGroupCollection *gc, dgConservationLaw *law);
+ static void registerBindings(binding *b);
};