From b21cb2e63fcbfacd91542325f9cdee2c14d6dbf8 Mon Sep 17 00:00:00 2001
From: Richard Comblen <richard.comblen@uclouvain.be>
Date: Thu, 4 Mar 2010 09:58:53 +0000
Subject: [PATCH] Tools to integrate a function
---
Common/LuaBindings.cpp | 2 +
Solver/CMakeLists.txt | 1 +
Solver/TESTCASES/MultirateAdvection.lua | 11 ++++
Solver/dgFunctionIntegrator.cpp | 53 +++++++++++++++++
Solver/dgFunctionIntegrator.h | 10 ++++
Solver/dgGroupOfElements.cpp | 78 +++++++++++++++++--------
6 files changed, 131 insertions(+), 24 deletions(-)
create mode 100644 Solver/dgFunctionIntegrator.cpp
create mode 100644 Solver/dgFunctionIntegrator.h
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index d5dd74b077..ba9f190a70 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -23,6 +23,7 @@
#include "dgRungeKuttaMultirate.h"
#include "dgSystemOfEquations.h"
#include "dgLimiter.h"
+#include "dgFunctionIntegrator.h"
#include "Bindings.h"
#include "dgResidual.h"
@@ -347,6 +348,7 @@ binding::binding(){
dgRungeKuttaMultirate43::registerBindings(this);
dgSlopeLimiterRegisterBindings(this);
dgSystemOfEquations::registerBindings(this);
+ dgFunctionIntegrator::registerBindings(this);
fullMatrix<double>::registerBindings(this);
function::registerBindings(this);
function::registerDefaultFunctions();
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index 1ce5fa3a7b..6aa81f0fe1 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -29,6 +29,7 @@ set(SRC
functionDerivator.cpp
luaFunction.cpp
functionSpace.cpp
+ dgFunctionIntegrator.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_SOURCE_DIR}/Solver *.h)
diff --git a/Solver/TESTCASES/MultirateAdvection.lua b/Solver/TESTCASES/MultirateAdvection.lua
index 6ec63c036e..47636aa8ba 100644
--- a/Solver/TESTCASES/MultirateAdvection.lua
+++ b/Solver/TESTCASES/MultirateAdvection.lua
@@ -20,6 +20,13 @@ function velocity( XYZ, FCT )
end
end
+function toIntegrate( XYZ, FCT )
+ for i=0,XYZ:size1()-1 do
+ X = XYZ:get(i,0)
+ Y = XYZ:get(i,1)
+ FCT:set(i,0,X*X)
+ end
+end
--[[
Example of a lua program driving the DG code
--]]
@@ -87,6 +94,7 @@ time=0
-- multirateRK:setLimiter(limiter)
--for i=1,1000
i=0
+integrator=dgFunctionIntegrator(functionLua(1, 'toIntegrate', {'XYZ'}):getName())
while time<0.2 do
-- norm1 = RK:iterate43SchlegelJCAM2009(law,dt,solution)
-- TEST with Explicit Euler multirate !!!
@@ -94,6 +102,9 @@ while time<0.2 do
time=time+dt
if (i % 10 == 0) then
print('*** ITER ***',i,time,norm2)
+ v=fullMatrix(1,1);
+ integrator:compute(solution2,v);
+ print('integral : ',v:get(0,0,0))
end
if (i % 10 == 0) then
-- solution:exportMsh(string.format("output/rt-%06d", i))
diff --git a/Solver/dgFunctionIntegrator.cpp b/Solver/dgFunctionIntegrator.cpp
new file mode 100644
index 0000000000..58822d6b87
--- /dev/null
+++ b/Solver/dgFunctionIntegrator.cpp
@@ -0,0 +1,53 @@
+#include "function.h"
+#include "dgFunctionIntegrator.h"
+#include "dgDofContainer.h"
+#include "fullMatrix.h"
+#include "dgGroupOfElements.h"
+
+#include <stdio.h>
+
+dgFunctionIntegrator::dgFunctionIntegrator(std::string fName):_fName(fName){}
+
+void dgFunctionIntegrator::compute(dgDofContainer *sol,fullMatrix<double> &result){
+ int nbFields=sol->getNbFields();
+ dataCacheMap cacheMap;
+ dataCacheDouble &UVW=cacheMap.provideData("UVW", 1, 3);
+ dataCacheDouble &solutionQPe=cacheMap.provideData("Solution", 1, nbFields);
+ dataCacheElement &cacheElement=cacheMap.getElement();
+ function *f=function::get(_fName);
+ dataCacheDouble *F=f->newDataCache(&cacheMap);
+ int nbRowResult=result.size1();
+ result.scale(0.0);
+ for(int iGroup=0;iGroup<sol->getGroups()->getNbElementGroups();iGroup++){
+ dgGroupOfElements &group=*sol->getGroups()->getElementGroup(iGroup);
+ fullMatrix<double> &solProxy=sol->getGroupProxy(&group);
+ cacheMap.setNbEvaluationPoints(group.getNbIntegrationPoints());
+ UVW.set(group.getIntegrationPointsMatrix());
+ fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
+ group.getCollocationMatrix().mult(solProxy , solutionQP);
+ fullMatrix<double> IPMatrix = group.getIntegrationPointsMatrix();
+ for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
+ solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
+ cacheElement.set(group.getElement(iElement));
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iElement, iPt);
+ for (int k=0;k<nbRowResult;k++){
+ result(0,k) += (*F)(iPt,k)*detJ*IPMatrix(iPt,3);
+ }
+ }
+ }
+ }
+}
+
+#include "Bindings.h"
+void dgFunctionIntegrator::registerBindings(binding *b){
+ classBinding *cb = b->addClass<dgFunctionIntegrator>("dgFunctionIntegrator");
+ cb->setDescription("Integrates a function, using the compute function");
+ methodBinding *mb;
+ mb = cb->setConstructor<dgFunctionIntegrator,std::string>();
+ mb->setArgNames("functionName",NULL);
+ mb->setDescription("a new dgFunctionIntegrator, get the solution using the compute method");
+ mb = cb->addMethod("compute", &dgFunctionIntegrator::compute);
+ mb->setArgNames("dofContainer","result",NULL);
+ mb->setDescription("compute the integral of the function");
+}
diff --git a/Solver/dgFunctionIntegrator.h b/Solver/dgFunctionIntegrator.h
new file mode 100644
index 0000000000..ebbb803764
--- /dev/null
+++ b/Solver/dgFunctionIntegrator.h
@@ -0,0 +1,10 @@
+#include <string>
+class dgDofContainer;
+class binding;
+class dgFunctionIntegrator{
+ std::string _fName;
+ public:
+ dgFunctionIntegrator(std::string fName);
+ void compute(dgDofContainer *sol,fullMatrix<double> &result);
+ static void registerBindings(binding *b);
+};
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index b9aceb8445..854e23e1aa 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1052,6 +1052,12 @@ void dgGroupCollection::buildGroupsOfInterfaces()
// Split the groups of elements depending on their local time step
double dgGroupCollection::splitGroupsForMultirate(int maxLevels,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
+ // outerBuffer: elements that use the big time step but talks to elements using the small time step
+
+
Msg::Info("Splitting Groups for multirate time stepping");
maxLevels--;// Number becomes maximum id
@@ -1059,8 +1065,14 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
std::vector<dgMiniInterface> *miniInterfaceV=_createMiniInterfaces(*this);
// elementToNeighbors[oldGroupId][oldElementId][neighborId (i.e. 0 to 2 for triangles)]<neighborOldGroupId,neighborOldElementId>
std::vector<std::vector<std::vector<std::pair<int,int> > > >elementToNeighbors;
-
+ // id of the new group, indexed by [oldGroupId][oldElementId]
+ // -1 if the element is not yet in a new group
+ // -2 if the element is in a new group whose id is not determined yet
std::vector<std::vector<int> >newGroupIds;
+ std::vector< dgGroupOfElements* >newGroups;// indexed by newGroupId
+ // localDt[oldGroupId][oldElementId]
+ std::vector<std::vector<double> >localDt;
+
newGroupIds.resize(getNbElementGroups());
elementToNeighbors.resize(getNbElementGroups());
for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
@@ -1068,6 +1080,8 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
elementToNeighbors[iGroup].resize(g->getNbElements());
newGroupIds[iGroup].assign(g->getNbElements(),-1);
}
+
+ // Build elementToNeighbors table (needed to have random access to the neighbors)
for(int iInterface=0;iInterface<miniInterfaceV->size();iInterface++){
dgMiniInterface &interface=miniInterfaceV->at(iInterface);
for(int iConn=0;iConn<interface.connectedElements.size();iConn++){
@@ -1082,10 +1096,9 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
}
}
- // find the range of time steps
+ // Compute the time step constrain per element
double dtMin = DBL_MAX;
double dtMax = 0;
- std::vector<std::vector<double> >localDt;
localDt.resize(getNbElementGroups());
for (int i=0;i<getNbElementGroups();i++){
dgAlgorithm::computeElementaryTimeSteps(*claw, *getElementGroup(i), solution->getGroupProxy(i), localDt[i]);
@@ -1102,12 +1115,16 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
MPI_Allreduce((void *)&dtMax, &dtMax_max, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
dtMax=dtMax_max;
#endif
- // dtMin is the time step for the most constrained element.
Msg::Info("Time step for standard RK should be %e",dtMin);
Msg::Info("Multirate base time step should be %e",dtMax);
+
+
+ // We take the base time step as 80% of the largest elementary time-step.
+ // This is arbitrary and should be optimized (maybe on the fly...)
double dtRef=dtMax*0.8;
+
// time steps are dtRef*2^(-dtExponent), with dtExponent ranging in [0:dtMaxExponent]
int dtMaxExponent=0;
while(dtRef/pow(2.0,(double)dtMaxExponent)>dtMin)
@@ -1117,30 +1134,38 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
dtMaxExponent=maxLevels;
}
_dtMaxExponent=dtMaxExponent;
- std::vector<MElement *>currentNewGroup;
- std::vector< dgGroupOfElements* >newGroups;// indexed by newGroupId
+
+ // Bounds to be able to access elements at the previous level (to add the neighbors in the current group)
int lowerLevelGroupIdStart=-1;
int lowerLevelGroupIdEnd=-1;
bool isOuterBufferLayer=false;
int currentNewGroupId=0;
int loopId=0;
+
+ // We pass two times by each exponent level: one for the normal(bulk) + innerBuffer groups, and one for the outerBuffer groups
for(int currentExponent=dtMaxExponent;currentExponent>=0;(!isOuterBufferLayer)?currentExponent--:currentExponent=currentExponent){
double currentDt=dtRef/pow(2.0,(double)currentExponent);
- std::map<int,std::vector<std::pair<int,int> > >mapNewGroups;
+ // element for the new groups at this level:
+ // mapNewGroups[oldGroupId][oldElementId]
+ std::map<int,std::vector<int> >mapNewGroups;
if(lowerLevelGroupIdStart==-1){
lowerLevelGroupIdStart=0;
}
else{
// Add the neighbors elements to the new groups
+ // For buffer AND non buffer layers
int _lowerLevelGroupIdStart=lowerLevelGroupIdStart;
int _lowerLevelGroupIdEnd=lowerLevelGroupIdEnd;
lowerLevelGroupIdStart=lowerLevelGroupIdEnd;
for(int iInterface=0;iInterface<miniInterfaceV->size();iInterface++){
dgMiniInterface &interface=miniInterfaceV->at(iInterface);
bool toAdd=false;
+ // if one of the elements adjacent to the interface is mapped to the previous level
+ // we check all elements adjacent to this interface, and add those that does not
+ // already have a new group to the current new group
for(int iConn=0;iConn<interface.connectedElements.size();iConn++){
int gId=interface.connectedElements[iConn].first;
int eId=interface.connectedElements[iConn].second;
@@ -1156,7 +1181,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
int eId=interface.connectedElements[iConn].second;
int newGroupId=newGroupIds[gId][eId];
if(newGroupId==-1){
- mapNewGroups[gId].push_back(std::pair<int,int>(gId,eId));
+ mapNewGroups[gId].push_back(eId);
newGroupIds[gId][eId]=-2;
}
}
@@ -1164,27 +1189,29 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
}
}
if(!isOuterBufferLayer){
+ // We add all the elements that are stable with the current time step, but not with twice that time step
for(int iGroup=0;iGroup<getNbElementGroups();iGroup++){
dgGroupOfElements *elGroup=getElementGroup(iGroup);
for(int iElement=0;iElement<elGroup->getNbElements();iElement++){
if(localDt[iGroup][iElement]>=currentDt && (localDt[iGroup][iElement]<currentDt*2 || currentExponent==0)){
if(newGroupIds[iGroup][iElement]==-1){
- mapNewGroups[iGroup].push_back(std::pair<int,int>(iGroup,iElement));
+ mapNewGroups[iGroup].push_back(iElement);
newGroupIds[iGroup][iElement]=-2;
}
}
}
}
+
+ // We split those elements into bulk and innerBuffer (i.e. those who have a neighbor with a bigger time step)
lowerLevelGroupIdStart=currentNewGroupId;
- for (std::map<int, std::vector<std::pair<int,int> > >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
+ for (std::map<int, std::vector<int> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
if(!it->second.empty()){
- std::vector<std::pair<int,int> >forBulk;
- std::vector<std::pair<int,int> >forInnerBuffer;
+ std::vector<int>forBulk;
+ std::vector<int>forInnerBuffer;
for(int i=0;i<it->second.size();i++){
bool inInnerBuffer=false;
- //std::vector<std::vector<std::vector<std::pair<int,int> > > >elementToNeighbors;
- int oldGId=it->second[i].first;
- int oldEId=it->second[i].second;
+ int oldGId=it->first;
+ int oldEId=it->second[i];
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){
@@ -1193,21 +1220,21 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
}
}
if(inInnerBuffer){
- forInnerBuffer.push_back(std::pair<int,int>(oldGId,oldEId));
+ forInnerBuffer.push_back(oldEId);
}
else{
- forBulk.push_back(std::pair<int,int>(oldGId,oldEId));
+ forBulk.push_back(oldEId);
}
}
for(int i=0;i<forBulk.size();i++){
- newGroupIds[forBulk[i].first][forBulk[i].second]=currentNewGroupId;
+ newGroupIds[it->first][forBulk[i]]=currentNewGroupId;
}
dgGroupOfElements *oldGroup=getElementGroup(it->first);
dgGroupOfElements *newGroup;
if(!forBulk.empty()){
std::vector<MElement*>forBulkV;
for(int i=0;i<forBulk.size();i++){
- forBulkV.push_back(getElementGroup(forBulk[i].first)->getElement(forBulk[i].second));
+ forBulkV.push_back(getElementGroup(it->first)->getElement(forBulk[i]));
}
newGroup=new dgGroupOfElements(forBulkV,oldGroup->getOrder(),oldGroup->getGhostPartition());
newGroup->copyPrivateDataFrom(oldGroup);
@@ -1222,7 +1249,7 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
if(!forInnerBuffer.empty()){
std::vector<MElement*>forInnerBufferV;
for(int i=0;i<forInnerBuffer.size();i++){
- forInnerBufferV.push_back(getElementGroup(forInnerBuffer[i].first)->getElement(forInnerBuffer[i].second));
+ forInnerBufferV.push_back(getElementGroup(it->first)->getElement(forInnerBuffer[i]));
}
newGroup=new dgGroupOfElements(forInnerBufferV,oldGroup->getOrder(),oldGroup->getGhostPartition());
newGroup->copyPrivateDataFrom(oldGroup);
@@ -1239,17 +1266,17 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
}
}
else{
- for (std::map<int, std::vector<std::pair<int,int> > >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
+ for (std::map<int, std::vector<int> >::iterator it = mapNewGroups.begin(); it !=mapNewGroups.end() ; ++it){
if(!it->second.empty()){
for(int i=0;i<it->second.size();i++){
- int oldGId=it->second[i].first;
- int oldEId=it->second[i].second;
+ int oldGId=it->first;
+ int oldEId=it->second[i];
newGroupIds[oldGId][oldEId]=currentNewGroupId;
}
dgGroupOfElements *oldGroup=getElementGroup(it->first);
std::vector<MElement *>newGroupV;
for(int i=0;i<it->second.size();i++){
- newGroupV.push_back(getElementGroup(it->second[i].first)->getElement(it->second[i].second));
+ newGroupV.push_back(getElementGroup(it->first)->getElement(it->second[i]));
}
dgGroupOfElements *newGroup=new dgGroupOfElements(newGroupV,oldGroup->getOrder(),oldGroup->getGhostPartition());
newGroup->copyPrivateDataFrom(oldGroup);
@@ -1267,6 +1294,9 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
lowerLevelGroupIdEnd=currentNewGroupId;
isOuterBufferLayer=!isOuterBufferLayer;
}
+
+
+ // Some stats
int count=0;
for(int i=0;i<newGroups.size();i++){
Msg::Info("%d New group # %d has %d elements",newGroups[i]->getMultirateExponent(),i,newGroups[i]->getNbElements());
--
GitLab