From 3008601ad2d09a0a2fc95635b2d226e0530c7f59 Mon Sep 17 00:00:00 2001
From: Bruno Seny <bruno.seny@student.uclouvain.be>
Date: Fri, 19 Mar 2010 14:01:36 +0000
Subject: [PATCH] Splitting of the groups for multirate performed on multirate
object
---
.../TESTCASES/ConservationTestMRAdvection.lua | 167 +++++++++++-------
Solver/dgRungeKuttaMultirate.cpp | 36 +++-
Solver/dgRungeKuttaMultirate.h | 46 ++---
3 files changed, 163 insertions(+), 86 deletions(-)
diff --git a/Solver/TESTCASES/ConservationTestMRAdvection.lua b/Solver/TESTCASES/ConservationTestMRAdvection.lua
index 734b937a0e..32cb24bd49 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 65e5142fe3..5c2392f396 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 100cac8203..9aadec32ca 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);
};
--
GitLab