Skip to content
Snippets Groups Projects
Commit f787456e authored by Richard Comblen's avatar Richard Comblen
Browse files

Multirate RK works with order 3 in all case, added possibility to define... 

Multirate RK works with order 3 in all case, added possibility to define buffer layer of several elements (but this is not needed yet)
parent 34f3d761
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 291 additions and 114 deletions
Solver/TESTCASES/ConvergenceMultirate.lua
+ 4
− 3
View file @ f787456e
......@@ -17,8 +17,8 @@ function velocity( XYZ, FCT )
for i=0,XYZ:size1()-1 do
X = XYZ:get(i,0)
Y = XYZ:get(i,1)
FCT:set(i,0,1)
FCT:set(i,1,0)
FCT:set(i,0,0)
FCT:set(i,1,1)
FCT:set(i,2,0)
end
end
......@@ -62,12 +62,13 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
GC=dgGroupCollection(myModel,2,order)
solTmp=dgDofContainer(GC,1)
solTmp:L2Projection(FS)
dt=GC:splitGroupsForMultirate(3,law,solTmp)
dt=GC:splitGroupsForMultirate(20,1,law,solTmp)
GC:buildGroupsOfInterfaces(myModel,2,order)
solution=dgDofContainer(GC,1)
solutionRef=dgDofContainer(GC,1)
solutionRef:L2Projection(FS)
solution:exportGroupIdMsh()
solution:exportMultirateGroupMsh()
solutionRef:exportMsh("solution")
nRefSteps=100
......
Solver/TESTCASES/MultirateAdvection.lua
+ 4
− 3
View file @ f787456e
......@@ -17,6 +17,7 @@ function velocity( XYZ, FCT )
Y = XYZ:get(i,1)
FCT:set(i,0,0)
FCT:set(i,1,1)
FCT:set(i,2,0)
end
end
......@@ -48,7 +49,7 @@ elseif (order == 5) then
end
print'*** Create a dg solver ***'
velF = functionLua(2, 'velocity', {'XYZ'}):getName()
velF = functionLua(3, 'velocity', {'XYZ'}):getName()
law=dgConservationLawAdvectionDiffusion(velF,"")
zero=functionConstant({0.}):getName();
......@@ -59,7 +60,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
GC=dgGroupCollection(myModel,2,order)
solTmp=dgDofContainer(GC,1)
solTmp:L2Projection(FS)
dt=GC:splitGroupsForMultirate(4,law,solTmp)
dt=GC:splitGroupsForMultirate(2,10,law,solTmp)
GC:buildGroupsOfInterfaces(myModel,2,order)
solution=dgDofContainer(GC,1)
solution:L2Projection(FS)
......@@ -95,7 +96,7 @@ time=0
--for i=1,1000
i=0
integrator=dgFunctionIntegrator(functionLua(1, 'toIntegrate', {'XYZ'}):getName())
while time<0.2 do
while i<1 do
-- norm1 = RK:iterate43SchlegelJCAM2009(law,dt,solution)
-- TEST with Explicit Euler multirate !!!
norm2 = multirateRK:iterate(dt,solution2)
......
Solver/dgDofContainer.cpp
+ 87
− 2
View file @ f787456e
......@@ -192,6 +192,19 @@ double dgDofContainer::norm() {
return localNorm;
#endif
}
double dgDofContainer::norm(std::vector<dgGroupOfElements *>gV){
double norm=0;
for(int i=0;i<gV.size();i++){
dgGroupOfElements* g=gV[i];
fullMatrix<double>&groupData=getGroupProxy(g);
for(int r=0;r<groupData.size1();r++){
for(int c=0;c<groupData.size1();c++){
norm+=groupData(r,c)*groupData(r,c);
}
}
}
return norm;
}
void dgDofContainer::save(const std::string name) {
FILE *f = fopen (name.c_str(),"wb");
_data->binarySave(f);
......@@ -233,7 +246,7 @@ void dgDofContainer::exportGroupIdMsh(){
// the elementnodedata::export does not work !!
std::ostringstream name_oss;
name_oss<<"groups.msh";
name_oss<<"groupsId.msh";
if(Msg::GetCommSize()>1)
name_oss<<"_"<<Msg::GetCommRank();
FILE *f = fopen (name_oss.str().c_str(),"w");
......@@ -294,6 +307,76 @@ void dgDofContainer::exportGroupIdMsh(){
*/
}
void dgDofContainer::exportMultirateGroupMsh(){
// the elementnodedata::export does not work !!
std::ostringstream name_oss;
name_oss<<"groupsMultirateType.msh";
if(Msg::GetCommSize()>1)
name_oss<<"_"<<Msg::GetCommRank();
FILE *f = fopen (name_oss.str().c_str(),"w");
int COUNT = 0;
for (int i=0;i < _groups.getNbElementGroups() ;i++){
COUNT += _groups.getElementGroup(i)->getNbElements();
}
fprintf(f,"$MeshFormat\n2.1 0 8\n$EndMeshFormat\n");
fprintf(f,"$ElementNodeData\n");
fprintf(f,"1\n");
fprintf(f,"\"%s\"\n","GroupsIds");
fprintf(f,"1\n");
fprintf(f,"0.0\n");
fprintf(f,"%d\n", Msg::GetCommSize() > 1 ? 4 : 3);
fprintf(f,"0\n 1\n %d\n",COUNT);
if(Msg::GetCommSize() > 1) fprintf(f,"%d\n", Msg::GetCommRank());
for (int i=0;i < _groups.getNbElementGroups() ;i++){
dgGroupOfElements *group = _groups.getElementGroup(i);
int groupType=0;//-3*group->getMultirateExponent();
if(group->getIsInnerMultirateBuffer())
groupType+=1;
if(group->getIsOuterMultirateBuffer())
groupType+=2;
for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
MElement *e = group->getElement(iElement);
int num = e->getNum();
fullMatrix<double> sol(getGroupProxy(i), iElement*_nbFields,_nbFields);
fprintf(f,"%d %d",num,sol.size1());
for (int k=0;k<sol.size1();++k) {
fprintf(f," %.16E ",groupType*1.0);
}
fprintf(f,"\n");
}
}
fprintf(f,"$EndElementNodeData\n");
fclose(f);
return;
// we should discuss that : we do a copy of the solution, this should
// be avoided !
/*std::map<int, std::vector<double> > data;
for (int i=0;i < _groups.getNbElementGroups() ;i++){
dgGroupOfElements *group = _groups.getElementGroup(i);
for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
MElement *e = group->getElement(iElement);
int num = e->getNum();
fullMatrix<double> sol(getGroupProxy(i), iElement*_nbFields,_nbFields);
std::vector<double> val;
// val.resize(sol.size2()*sol.size1());
val.resize(sol.size1());
int counter = 0;
// for (int iC=0;iC<sol.size2();iC++)
printf("%g %g %g\n",sol(0,0),sol(1,0),sol(2,0));
for (int iR=0;iR<sol.size1();iR++)val[counter++] = sol(iR,0);
data[num] = val;
}
}
PView *pv = new PView (name, "ElementNodeData", _gm, data, 0.0, 1);
pv->getData()->writeMSH(name+".msh", false);
delete pv;
*/
}
void dgDofContainer::exportMsh(const std::string name)
{
// the elementnodedata::export does not work !!
......@@ -386,7 +469,9 @@ void dgDofContainer::registerBindings(binding *b){
cm->setArgNames("filename", NULL);
cm = cb->addMethod("exportGroupIdMsh",&dgDofContainer::exportGroupIdMsh);
cm->setDescription("Export the group ids for gmsh visualization");
cm = cb->addMethod("norm",&dgDofContainer::norm);
cm = cb->addMethod("exportMultirateGroupMsh",&dgDofContainer::exportMultirateGroupMsh);
cm->setDescription("Export the group Multirate properties for gmsh visualization");
cm = cb->addMethod("norm",(double (dgDofContainer::*)())&dgDofContainer::norm);
cm->setDescription("Returns the norm of the vector");
cm = cb->addMethod("scale",(void (dgDofContainer::*)(double))&dgDofContainer::scale);
cm->setArgNames("factor",NULL);
......
Solver/dgDofContainer.h
+ 2
− 0
View file @ f787456e
......@@ -28,6 +28,7 @@ public:
void scale(double f);
void scale(std::vector<dgGroupOfElements*>groups, double f);
double norm();
double norm(std::vector<dgGroupOfElements*>groups);
void axpy(dgDofContainer &x, double a=1.);
void axpy(std::vector<dgGroupOfElements*>groups,dgDofContainer &x, double a=1.);
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
......@@ -44,6 +45,7 @@ public:
void Mesh2Mesh_L2Projection(dgDofContainer &other);
void exportMsh(const std::string name);
void exportGroupIdMsh();
void exportMultirateGroupMsh();
static void registerBindings(binding *b);
inline dgGroupCollection *getGroups() { return &_groups; }
......
Solver/dgGroupOfElements.cpp
+ 54
− 33
View file @ f787456e
......@@ -19,6 +19,7 @@
#include <string.h>
#endif
static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
int npts;
IntPt *pts;
......@@ -737,7 +738,7 @@ void dgGroupCollection::buildParallelStructure()
}
// Split the groups of elements depending on their local time step
double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution) {
double dgGroupCollection::splitGroupsForMultirate(int maxLevels,int bufferSize,dgConservationLaw *claw, dgDofContainer *solution) {
// What are the levels/layers:
// bulk: elements that are time stepped using the "normal" 4 stage Runge-Kutta
// innerBuffer: elements that use the small time step but talks to elements using the big time step
......@@ -845,9 +846,9 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
else{
// Add the neighbors elements to the new groups
// For buffer AND non buffer layers
int _lowerLevelGroupIdStart=lowerLevelGroupIdStart;
int _lowerLevelGroupIdEnd=lowerLevelGroupIdEnd;
lowerLevelGroupIdStart=lowerLevelGroupIdEnd;
// We add bufferSize elements (most of the time, bufferSize=1 is enough)
for(int iLoop=0;iLoop<bufferSize;iLoop++){
for(int iInterface=0;iInterface<miniInterfaceV->size();iInterface++){
dgMiniInterface &interface=miniInterfaceV->at(iInterface);
bool toAdd=false;
......@@ -858,7 +859,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
int gId=interface.connections[iConn].iGroup;
int eId=interface.connections[iConn].iElement;
int newGroupId=newGroupIds[gId][eId];
if(newGroupId>=0 /*newGroupId >= _lowerLevelGroupIdStart && newGroupId<_lowerLevelGroupIdEnd*/){
if(newGroupId>=0 || ( newGroupId>-2-iLoop && newGroupId<-1) ){
toAdd=true;
continue;
}
......@@ -870,7 +871,8 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
int newGroupId=newGroupIds[gId][eId];
if(newGroupId==-1){
mapNewGroups[gId].push_back(eId);
newGroupIds[gId][eId]=-2;
newGroupIds[gId][eId]=-2-iLoop;
}
}
}
}
......@@ -892,22 +894,41 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
// We split those elements into bulk and innerBuffer (i.e. those who have a neighbor with a bigger time step)
lowerLevelGroupIdStart=currentNewGroupId;
#define MAXBUFFERSIZE 100000
if(currentExponent>0){
for(int iLoop=0;iLoop<bufferSize;iLoop++){
for (std::map<int, std::vector<int> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
if(!it->second.empty()){
std::vector<int>forBulk;
std::vector<int>forInnerBuffer;
for(int i=0;i<it->second.size();i++){
bool inInnerBuffer=false;
int oldGId=it->first;
int oldEId=it->second[i];
// We only consider elements in the map and not in the inner buffer for a different iLoop
if(newGroupIds[oldGId][oldEId]>-2 || newGroupIds[oldGId][oldEId]<-(MAXBUFFERSIZE-1) )
continue;
for(int iNeighbor=0;iNeighbor<elementToNeighbors[oldGId][oldEId].size();iNeighbor++){
std::pair<int,int> neighIds=elementToNeighbors[oldGId][oldEId][iNeighbor];
if(newGroupIds[neighIds.first][neighIds.second]==-1){
if(newGroupIds[neighIds.first][neighIds.second]==-1-iLoop*MAXBUFFERSIZE){
inInnerBuffer=true;
continue;
}
}
if(inInnerBuffer){
newGroupIds[oldGId][oldEId]=-1-(iLoop+1)*MAXBUFFERSIZE;
}
}
}
}
}
for (std::map<int, std::vector<int> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
if(!it->second.empty()){
std::vector<int>forBulk;
std::vector<int>forInnerBuffer;
for(int i=0;i<it->second.size();i++){
int oldGId=it->first;
int oldEId=it->second[i];
if( newGroupIds[oldGId][oldEId]>-2 || newGroupIds[oldGId][oldEId]<-(MAXBUFFERSIZE-1) ){
forInnerBuffer.push_back(oldEId);
}
else{
......@@ -994,7 +1015,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
// Some stats
int count=0;
for(int i=0;i<newGroups.size();i++){
Msg::Info("old: %d, level %d, New group # %d has %d elements",oldGroupIds[i],newGroups[i]->getMultirateExponent(),i,newGroups[i]->getNbElements());
Msg::Info("Old group #%d, exponent %d, New group #%d has %d elements",oldGroupIds[i],newGroups[i]->getMultirateExponent(),i,newGroups[i]->getNbElements());
if(newGroups[i]->getIsInnerMultirateBuffer())
Msg::Info("InnerBuffer");
else if(newGroups[i]->getIsOuterMultirateBuffer())
......@@ -1182,7 +1203,7 @@ void dgGroupCollection::registerBindings(binding *b)
cm->setDescription("Build the group of interfaces, i.e. boundary interfaces and inter-element interfaces");
cm = cb->addMethod("splitGroupsForMultirate",&dgGroupCollection::splitGroupsForMultirate);
cm->setDescription("Split the groups according to their own stable time step");
cm->setArgNames("maxLevels","claw","solution",NULL);
cm->setArgNames("maxLevels","bufferSize","claw","solution",NULL);
cm = cb->addMethod("splitGroupsByVerticalLayer",&dgGroupCollection::splitGroupsByVerticalLayer);
cm->setDescription("Split the groups according vertical layer structure. The first is defined by the topLevelTags.");
cm->setArgNames("topLevelTags",NULL);
......
Solver/dgGroupOfElements.h
+ 1
− 1
View file @ f787456e
......@@ -247,7 +247,7 @@ class dgGroupCollection {
void buildGroupsOfElements (GModel *model,int dimension, int order);
void buildGroupsOfInterfaces ();
double splitGroupsForMultirate(int maxLevels,dgConservationLaw *claw, dgDofContainer *solution);
double splitGroupsForMultirate(int maxLevels,int bufferSize,dgConservationLaw *claw, dgDofContainer *solution);
void splitGroupsByVerticalLayer(std::vector<std::string> topLevelTags);
void find (MElement *elementToFind, int &iGroup, int &ithElementOfGroup);
......
Solver/dgRungeKuttaMultirate.cpp
+ 138
− 72
View file @ f787456e
......@@ -6,12 +6,56 @@
#include "dgGroupOfElements.h"
#include <algorithm>
#include <limits.h>
#include <stdio.h>
//#define MULTIRATEVERBOSE
static double A43[4][4]={
{0, 0, 0, 0},
{1.0/2.0, 0, 0 ,0},
{-1.0/6.0, 2.0/3.0, 0, 0},
{1.0/3.0, -1.0/3.0, 1, 0}
};
static double AInner43[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/4.0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, -1.0/12.0, 1.0/3.0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
};
// Big step RK43
static double AOuter43[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
};
static double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
static double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
static double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
static double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
static double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
static double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservationLaw *law, int nStages){
_law=law;
_residualVolume=new dgResidualVolume(*_law);
_residualInterface=new dgResidualInterface(*_law);
_residualBoundary=new dgResidualBoundary(*_law);
_gc=gc;
_K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
_K[i]=new dgDofContainer(gc,_law->getNbFields());
......@@ -40,12 +84,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
_bulkGroupsOfElements[ge->getMultirateExponent()].first.push_back(ge);
}
}
// std::vector<std::set<dgGroupOfFaces*> >bulkFacesSet;
// std::vector<std::set<dgGroupOfFaces*> >innerBufferFacesSet;
// std::vector<std::set<dgGroupOfFaces*> >outerBufferFacesSet;
// bulkFacesSet.resize(_maxExponent+1);
// innerBufferFacesSet.resize(_maxExponent+1);
// outerBufferFacesSet.resize(_maxExponent+1);
for(int iGroup=0;iGroup<gc->getNbFaceGroups();iGroup++){
dgGroupOfFaces *gf=gc->getFaceGroup(iGroup);
const dgGroupOfElements *ge[2];
......@@ -66,7 +104,6 @@ dgRungeKuttaMultirate::dgRungeKuttaMultirate(dgGroupCollection* gc,dgConservatio
dgGroupOfFaces *gf=gc->getBoundaryGroup(iGroup);
const dgGroupOfElements *ge[1];
ge[0]=&gf->getGroupOfConnections(0).getGroupOfElements();
//ge[1]=&gf->getGroupRight();
for(int i=0;i<1;i++){
if(ge[i]->getIsInnerMultirateBuffer()){
_innerBufferGroupsOfElements[ge[i]->getMultirateExponent()].second.push_back(gf);
......@@ -103,7 +140,11 @@ dgRungeKuttaMultirate::~dgRungeKuttaMultirate(){
}
void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
//Msg::Info("Compute K%d at level %d %s",iK,exponent,isBuffer?"Buffer":"Bulk");
#ifdef MULTIRATEVERBOSE
for(int i=0;i<exponent;i++)
printf("\t");
printf("Exponent %d, %s, compute K%d\n",exponent,isBuffer?" Buffer":"Not buffer",iK);
#endif
if(isBuffer){
std::vector<dgGroupOfElements *>&vei=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&veo=_outerBufferGroupsOfElements[exponent].first;
......@@ -124,7 +165,6 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
dgGroupOfFaces *faces=*it;
if(faces->getNbGroupOfConnections()==1){
_residualBoundary->computeAndMap1Group(*faces,*_currentInput,*_residual);
//Msg::Info("Buffer face group %p is boundary in multirate",faces);
}
else{
const dgGroupOfElements *gL = &faces->getGroupOfConnections(0).getGroupOfElements();
......@@ -133,13 +173,14 @@ void dgRungeKuttaMultirate::computeK(int iK,int exponent,bool isBuffer){
fullMatrix<double> residuInterface(faces->getNbNodes(),faces->getNbElements()*2*_law->getNbFields());
faces->mapToInterface(_law->getNbFields(), _currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR), solInterface);
_residualInterface->compute1Group(*faces,solInterface,_currentInput->getGroupProxy(gL), _currentInput->getGroupProxy(gR),residuInterface);
//Msg::Info("Buffer face group %p is mapped left or right in multirate",faces);
if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer())
if(gL->getMultirateExponent()==exponent && gL->getIsMultirateBuffer()){
faces->mapLeftFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gL));
if(gR->getMultirateExponent()==exponent && gR->getIsMultirateBuffer())
}
if(gR->getMultirateExponent()==exponent && gR->getIsMultirateBuffer()){
faces->mapRightFromInterface(_law->getNbFields(), residuInterface,_residual->getGroupProxy(gR));
}
}
}
fullMatrix<double> iMassEl, KEl,residuEl;
for(int iGroup=0;iGroup<vei.size();iGroup++){
dgGroupOfElements *group=vei[iGroup];
......@@ -236,43 +277,25 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
if(exponent>_maxExponent){
return;
}
//Msg::Info("Input for K%d at level %d %s",iK,exponent,isBuffer?"Buffer":"Bulk");
#ifdef MULTIRATEVERBOSE
for(int i=0;i<exponent;i++)
printf("\t");
printf("Exponent %d, %s, input K%d\n",exponent,isBuffer?" Buffer":"Not buffer",iK);
#endif
double localDt=_dt/pow(2.0,(double)exponent);
double _A[4][4]={
{0, 0, 0, 0},
{1.0/2.0, 0, 0 ,0},
{-1.0/6.0, 2.0/3.0, 0, 0},
{1.0/3.0, -1.0/3.0, 1, 0}
};
double _AInner[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 ,0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/12.0, 1.0/3.0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/4.0, 0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, -1.0/12.0, 1.0/3.0, 0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/6.0, -1.0/6.0, 1.0/2.0, 0, 0},
{1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0, 1.0/12.0, 1.0/6.0, 1.0/6.0, 1.0/12.0, 0}
};
// Big step RK43
double _AOuter[10][10]={
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/4.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1.0/2.0 , 0, 0, 0, 0, 0, 0, 0, 0, 0},
{-1.0/6.0, 0, 0, 0, 2.0/3.0, 0, 0, 0, 0, 0},
{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
{1.0/12.0, 0, 0, 0, 1.0/6.0, 1.0/2.0, 0, 0, 0, 0},
{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
{1.0/3.0 , 0, 0, 0, -1.0/3.0, 1, 0, 0, 0, 0},
};
// Msg::Info("K%d, %d, %s",iK,exponent,isBuffer?"Buffer":"Bulk");
// compute K[iK] for this exponent and buffer
if(isBuffer){
if(!isBuffer){
std::vector<dgGroupOfElements *> &gV=_bulkGroupsOfElements[exponent].first;
_currentInput->scale(gV,0.0);
_currentInput->axpy(gV,*_solution,1.0);
for(int i=0;i<iK;i++){
if(A43[iK][i]!=0.0){
_currentInput->axpy(gV,*_K[i],A43[iK][i]*localDt);
}
}
}
else{
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
_currentInput->scale(gVi,0.0);
......@@ -280,22 +303,45 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
_currentInput->axpy(gVi,*_solution,1.0);
_currentInput->axpy(gVo,*_solution,1.0);
for(int i=0;i<iK;i++){
if(_AInner[iK][i]!=0.0)
_currentInput->axpy(gVi,*_K[i],_AInner[iK][i]*localDt*2);
if(_AOuter[iK][i]!=0.0)
_currentInput->axpy(gVo,*_K[i],_AOuter[iK][i]*localDt*2);
if(AInner43[iK][i]!=0.0)
_currentInput->axpy(gVi,*_K[i],AInner43[iK][i]*localDt*2);
if(AOuter43[iK][i]!=0.0)
_currentInput->axpy(gVo,*_K[i],AOuter43[iK][i]*localDt*2);
}
/*
NOT NEEDED
// We need to update input for the neighboring elements with bigger time step.
// if there is no corresponding K to be computed
if( (iK>0 && iK<4) || (iK>5 && iK<9) ){
std::vector<dgGroupOfElements *>&gVbigger=_bulkGroupsOfElements[exponent-1].first;
_currentInput->scale(gVbigger,0.0);
_currentInput->axpy(gVbigger,*_solution,1.0);
int idOuter2Bulk[10]={0,0,0,0,1,2,2,2,2,3};
for(int i=0;i<iK;i++){
if(AOuter43[iK][i]!=0.0){
// ON DIRAIT QUE CETTE LIGNE DE CODE NE FAIT RIEN ...
_currentInput->axpy(gVbigger,*_K[idOuter2Bulk[i]],AOuter43[iK][i]*localDt*2);
}
}
else{
std::vector<dgGroupOfElements *> &gV=_bulkGroupsOfElements[exponent].first;
_currentInput->scale(gV,0.0);
_currentInput->axpy(gV,*_solution,1.0);
std::vector<dgGroupOfElements *>&gViBigger=_innerBufferGroupsOfElements[exponent-1].first;
_currentInput->scale(gViBigger,0.0);
_currentInput->axpy(gViBigger,*_solution,1.0);
// shift
int s=0;
if(upperLeveliK>=5){
for(int i=0;i<4;i++){
_currentInput->axpy(gViBigger,*_K[i],b[i]*localDt*2);
}
s+=5;
}
for(int i=0;i<iK;i++){
if(_A[iK][i]!=0.0){
_currentInput->axpy(gV,*_K[i],_A[iK][i]*localDt);
if(AOuter43[iK][i]!=0.0){
_currentInput->axpy(gViBigger,*_K[idOuter2Bulk[i]+s],AOuter43[iK][i]*localDt*2);
}
}
}
*/
}
if(!isBuffer){
switch(iK){
......@@ -318,12 +364,10 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
computeInputForK(iK%5,exponent,false);
}
}
computeK(iK,exponent,isBuffer);
// Msg::Info("Multirate %d %0.16e",iK,_K[iK]->norm());
if( (iK==3 && !isBuffer) || (iK==9 && isBuffer) ){
updateSolution(exponent,isBuffer);
}
if(!isBuffer){
if(exponent==0){
computeK(iK,exponent,false);
if(iK==3)
updateSolution(exponent,false);
switch(iK){
case 0:
for(int i=1;i<4;i++){
......@@ -337,17 +381,37 @@ void dgRungeKuttaMultirate43::computeInputForK(int iK,int exponent,bool isBuffer
break;
}
}
if(isBuffer && exponent>0){
if( (iK%5)<4)
computeK(iK%5,exponent,false);
computeK(iK,exponent,true);
if( (iK%5)==3)
updateSolution(exponent,false);
if(iK==9)
updateSolution(exponent,true);
switch(iK%5){
case 0:
for(int i=1;i<4;i++){
computeInputForK(i,exponent+1,true);
}
break;
case 2:
for(int i=6;i<9;i++){
computeInputForK(i,exponent+1,true);
}
break;
}
}
}
void dgRungeKuttaMultirate43::updateSolution(int exponent,bool isBuffer){
//Msg::Info("Updating solution at level %d %s",exponent,isBuffer?"Buffer":"Bulk");
#ifdef MULTIRATEVERBOSE
for(int i=0;i<exponent;i++)
printf("\t");
printf("Updating solution at level %d %s\n",exponent,isBuffer?"Buffer":"Bulk");
#endif
double localDt=_dt/pow(2.0,(double)exponent);
double b[4]={1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
double c[4]={0, 1.0/2.0, 1.0/2.0, 1};
double bInner[10]={1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0,1.0/12.0, 1.0/6.0, 1.0/6.0,1.0/12.0,0};
double cInner[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
double bOuter[10]={1.0/6.0, 0 , 0 , 0 , 1.0/3.0,1.0/3.0, 0 , 0 , 0 , 1.0/6.0};
double cOuter[10]={0, 1.0/4.0, 1.0/4.0, 1.0/2.0, 1.0/2.0, 1.0/2.0, 3.0/4.0, 3.0/4.0, 1, 1 };
if(isBuffer){
std::vector<dgGroupOfElements *>&gVi=_innerBufferGroupsOfElements[exponent].first;
std::vector<dgGroupOfElements *>&gVo=_outerBufferGroupsOfElements[exponent].first;
......@@ -364,6 +428,8 @@ void dgRungeKuttaMultirate43::updateSolution(int exponent,bool isBuffer){
if(b[i]!=0.0)
_solution->axpy(gV,*_K[i],b[i]*localDt);
}
_currentInput->scale(gV,0.0);
_currentInput->axpy(gV,*_solution,1.0);
}
}
......
Solver/dgRungeKuttaMultirate.h
+ 1
− 0
View file @ f787456e
......@@ -13,6 +13,7 @@ class dgRungeKuttaMultirate: public dgRungeKutta{
dgResidualVolume *_residualVolume;
dgResidualInterface *_residualInterface;
dgResidualBoundary *_residualBoundary;
dgGroupCollection *_gc;
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_bulkGroupsOfElements;// int is the multirateExponent
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_innerBufferGroupsOfElements;// int is the multirateExponent
std::vector<std::pair<std::vector<dgGroupOfElements*>,std::vector<dgGroupOfFaces*> > >_outerBufferGroupsOfElements;// int is the multirateExponent
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment