Skip to content
Snippets Groups Projects
Commit 3008601a authored by Bruno Seny's avatar Bruno Seny
Browse files

Splitting of the groups for multirate performed on multirate object

parent 7d74dd25
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Solver/TESTCASES/ConservationTestMRAdvection.lua
+ 107
− 60
View file @ 3008601a
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,26 +172,33 @@ 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))
......@@ -168,15 +206,24 @@ integrator=dgFunctionIntegrator(functionLua(1, 'toIntegrate', {functionSolution.
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)
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)
......
Solver/dgRungeKuttaMultirate.cpp
+ 31
− 5
View file @ 3008601a
......@@ -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,6 +960,9 @@ 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");
......
Solver/dgRungeKuttaMultirate.h
+ 25
− 21
View file @ 3008601a
......@@ -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,6 +43,7 @@ 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);
......@@ -60,6 +63,7 @@ class dgRungeKuttaMultirateConservative: 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);
dgRungeKuttaMultirateConservative(dgGroupCollection *gc,dgConservationLaw *law, fullMatrix<double> *A, fullMatrix<double> *b,fullMatrix<double> *c);
void printButcher();
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment