From c4529a00865bec678daf0fa79caae1a6729c13ed Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Tue, 16 Feb 2010 12:54:04 +0000
Subject: [PATCH] dg: cleaning structure of dgGroupOfElements and dgRungeKutta,
dgSystemOfEquations can mostly be removed (only limiter left)
---
Common/LuaBindings.cpp | 10 +-
Common/LuaBindings.h | 1 +
Solver/CMakeLists.txt | 1 +
Solver/TESTCASES/CylinderEddies.lua | 4 +-
Solver/TESTCASES/RayleighTaylor.lua | 18 +-
Solver/TESTCASES/Stommel.lua | 20 +-
Solver/TESTCASES/rect.geo | 4 +-
Solver/TESTCASES/supra.lua | 2 +
Solver/dgAlgorithm.cpp | 85 ++-----
Solver/dgAlgorithm.h | 25 +-
Solver/dgConservationLaw.cpp | 14 +-
Solver/dgConservationLaw.h | 2 +-
Solver/dgDofContainer.cpp | 122 +++++++++-
Solver/dgDofContainer.h | 7 +-
Solver/dgGroupOfElements.cpp | 359 ++++++++--------------------
Solver/dgGroupOfElements.h | 17 +-
Solver/dgRungeKutta.cpp | 90 +++++++
Solver/dgRungeKutta.h | 19 ++
Solver/dgSlopeLimiter.cpp | 7 +-
Solver/dgSystemOfEquations.cpp | 141 ++---------
Solver/dgSystemOfEquations.h | 18 +-
contrib/kbipack/gmp_matrix_io.c | 2 +-
22 files changed, 444 insertions(+), 524 deletions(-)
create mode 100644 Solver/dgRungeKutta.cpp
create mode 100644 Solver/dgRungeKutta.h
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index 1996015a45..71157cadbc 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -9,14 +9,17 @@
#include "GFace.h"
#include "DivideAndConquer.h"
#include "Bindings.h"
-#include "dgSystemOfEquations.h"
#include "luaFunction.h"
#include "function.h"
+#include "GModel.h"
#include "dgGroupOfElements.h"
+#include "dgDofContainer.h"
#include "dgConservationLawShallowWater2d.h"
#include "dgConservationLawAdvection.h"
#include "dgConservationLawPerfectGas.h"
#include "dgConservationLawWaveEquation.h"
+#include "dgRungeKutta.h"
+#include "dgSystemOfEquations.h"
extern "C" {
#include "lua.h"
@@ -262,13 +265,16 @@ binding::binding(){
GModel::registerBindings(this);
fullMatrix<double>::registerBindings(this);
function::registerBindings(this);
+ dgGroupCollection::registerBindings(this);
+ dgDofContainer::registerBindings(this);
dgConservationLaw::registerBindings(this);
- dgSystemOfEquations::registerBindings(this);
+ dgRungeKutta::registerBindings(this);
dgBoundaryCondition::registerBindings(this);
dgConservationLawShallowWater2dRegisterBindings(this);
dgConservationLawWaveEquationRegisterBindings(this);
dgConservationLawAdvectionDiffusionRegisterBindings(this);
dgPerfectGasLaw2dRegisterBindings(this);
+ dgSystemOfEquations::registerBindings(this);
functionLua::registerBindings(this);
function::registerDefaultFunctions();
MVertex::registerBindings(this);
diff --git a/Common/LuaBindings.h b/Common/LuaBindings.h
index c94c75f3b1..866a5058ca 100644
--- a/Common/LuaBindings.h
+++ b/Common/LuaBindings.h
@@ -27,6 +27,7 @@ class className{
}
static const std::string &get(){
if(_name==""){
+ Msg::Error("Class unknown to Lua");
throw;
}
return _name;
diff --git a/Solver/CMakeLists.txt b/Solver/CMakeLists.txt
index e6a12f90cf..80c9b312e0 100644
--- a/Solver/CMakeLists.txt
+++ b/Solver/CMakeLists.txt
@@ -13,6 +13,7 @@ set(SRC
multiscaleLaplace.cpp
dgGroupOfElements.cpp
dgAlgorithm.cpp
+ dgRungeKutta.cpp
dgConservationLaw.cpp
dgConservationLawAdvection.cpp
dgSystemOfEquations.cpp
diff --git a/Solver/TESTCASES/CylinderEddies.lua b/Solver/TESTCASES/CylinderEddies.lua
index b35dc895fe..33c66654ef 100644
--- a/Solver/TESTCASES/CylinderEddies.lua
+++ b/Solver/TESTCASES/CylinderEddies.lua
@@ -58,11 +58,11 @@ DG:exportSolution('output/cyl_0')
print'*** solve ***'
-CFL = 20.1;
+CFL = 5.1;
dt = CFL * DG:computeInvSpectralRadius();
print('DT = ',dt)
T = 0;
-for i=1,2 do
+for i=1,100000 do
dt = CFL * DG:computeInvSpectralRadius();
norm = DG:RK44(dt)
T = T + dt
diff --git a/Solver/TESTCASES/RayleighTaylor.lua b/Solver/TESTCASES/RayleighTaylor.lua
index dfd59a504b..4fa27208d5 100644
--- a/Solver/TESTCASES/RayleighTaylor.lua
+++ b/Solver/TESTCASES/RayleighTaylor.lua
@@ -48,14 +48,14 @@ if (order == 1) then
elseif (order == 2) then
myModel:load ('rect2.msh')
elseif (order == 3) then
- myModel:load ('rect2.msh')
+ myModel:load ('rect3.msh')
elseif (order == 4) then
myModel:load ('rect4.msh')
+elseif (order == 5) then
+ myModel:load ('rect5.msh')
end
print'*** Create a dg solver ***'
-DG = dgSystemOfEquations (myModel)
-DG:setOrder(order)
law=dgPerfectGasLaw2d()
g=fullMatrix(4,1);
@@ -65,18 +65,18 @@ g:set(2,0,-1.)
g:set(3,0,0)
law:setSource(functionConstant(g):getName())
-DG:setConservationLaw(law)
law:addBoundaryCondition('Walls',law:newSlipWallBoundary())
FS = functionLua(4, 'initial_condition', {'XYZ'}):getName()
law:addBoundaryCondition('Top',law:newOutsideValueBoundary(FS))
-DG:setup()
-
-print'*** setting the initial solution ***'
-
-DG:L2Projection(FS)
+GC=dgGroupCollection(myModel,2,order)
+solution=dgDofContainer(GC,4)
+solution:L2Projection(FS)
DG:limitSolution()
+GC:splitGroupsForMultirate(0.0003,law)
+GC:buildGroupsOfInterfaces(myModel,2,order)
+print'*** setting the initial solution ***'
--print'*** export ***'
DG:exportSolution('output/rt_0')
diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua
index 4342a7e812..e11f1cd760 100644
--- a/Solver/TESTCASES/Stommel.lua
+++ b/Solver/TESTCASES/Stommel.lua
@@ -1,14 +1,20 @@
model = GModel()
model:load ('stommel_square.msh')
-dg = dgSystemOfEquations (model)
order=3
-dg:setOrder (order)
+dimension = 2
claw = dgConservationLawShallowWater2d()
claw:addBoundaryCondition('Wall',claw:newBoundaryWall())
-dg:setConservationLaw(claw)
-dg:setup()
-
-dg:exportSolution('output/init')
+groups = dgGroupCollection(model, dimension, order)
+-- groups:split...
+groups:buildGroupsOfInterfaces()
+solution = dgDofContainer(groups, claw:getNbFields())
+a = fullMatrix(3,1);
+a:set(0,0, 1.);
+a:set(1,0, 2.);
+a:set(2,0, 3.);
+f = functionConstant(a):getName()
+solution:L2Projection(f)
+solution:exportMsh('output/init')
for i=1,100000 do
norm = dg:RK44(150*(3/(2.*order+1)))
@@ -16,6 +22,6 @@ for i=1,100000 do
print ('iter ', i, norm)
end
if ( i%100 ==0 ) then
- dg:exportSolution(string.format('output/solution-%06d',i))
+ solution:exportMsh(string.format('output/solution-%06d',i))
end
end
diff --git a/Solver/TESTCASES/rect.geo b/Solver/TESTCASES/rect.geo
index 7deb863a25..9ab82d5364 100644
--- a/Solver/TESTCASES/rect.geo
+++ b/Solver/TESTCASES/rect.geo
@@ -18,8 +18,8 @@ Plane Surface(11) = {10};
Physical Surface("sprut") = {11, 9};
Physical Line("Walls") = {5, 7, 3, 4, 1};
Physical Line("Top") = {6};
-Transfinite Line {5, 7, 1, 3} = 100 Using Progression 1;
-Transfinite Line {6, 4, 2} = 50 Using Progression 1;
+Transfinite Line {5, 7, 1, 3} = 50 Using Progression 1;
+Transfinite Line {6, 4, 2} = 25 Using Progression 1;
Transfinite Surface {9} Alternated;
Transfinite Surface {11} Alternated;
Recombine Surface {9, 11};
diff --git a/Solver/TESTCASES/supra.lua b/Solver/TESTCASES/supra.lua
index 7ad6ccc051..9187b68ad3 100644
--- a/Solver/TESTCASES/supra.lua
+++ b/Solver/TESTCASES/supra.lua
@@ -29,7 +29,9 @@ function valueLeft(f)
end
law = dgConservationLawAdvectionDiffusion('',functionLua(1,'diffusivity',{'Solution'}):getName())
+print(law)
dg:setConservationLaw(law)
+print(law)
-- boundary condition
--[[
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index db458c85e6..99fc7bce66 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -25,7 +25,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
{
// ----- 1 ---- get the solution at quadrature points
// ----- 1.1 --- allocate a matrix of size (nbFields * nbElements, nbQuadraturePoints)
- int nbFields = claw.nbFields();
+ int nbFields = claw.getNbFields();
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
// ----- 1.2 --- multiply the solution by the collocation matrix
group.getCollocationMatrix().mult(solution , solutionQP);
@@ -96,7 +96,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
}
if (sourceTerm){
- source.setAsProxy(Source, iElement*claw.nbFields(),claw.nbFields());
+ source.setAsProxy(Source, iElement*claw.getNbFields(),claw.getNbFields());
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iElement, iPt);
for (int k=0;k<nbFields;k++){
@@ -127,7 +127,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
{
// A) global operations before entering the loop over the faces
// A1 ) copy locally some references from the group for the sake of lisibility
- int nbFields = claw.nbFields();
+ int nbFields = claw.getNbFields();
int nbNodesLeft = group.getGroupLeft().getNbNodes();
int nbNodesRight = group.getGroupRight().getNbNodes();
int nbFaces = group.getNbElements();
@@ -254,66 +254,8 @@ void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
}
-void dgAlgorithm::rungeKutta (const dgConservationLaw &claw, // conservation law
- dgGroupCollection &groups,
- double h, // time-step
- dgDofContainer &sol,
- dgDofContainer &resd,
- dgLimiter *limiter,
- int orderRK
- ) // order of RK integrator
-{
-
- // U_{n+1}=U_n+h*(SUM_i a_i*K_i)
- // K_i=M^(-1)R(U_n+b_i*K_{i-1})
-
- double a[4] = {h/6.0,h/3.0,h/3.0,h/6.0};
- double b[4] = {0.,h/2.0,h/2.0,h};
-
- if (orderRK == 1){
- a[0] = h;
- }
-
- // Current updated solution
- fullMatrix<double> residuEl, KEl;
- fullMatrix<double> iMassEl;
-
- int nbFields = claw.nbFields();
- dgDofContainer K (groups,nbFields);
- dgDofContainer Unp (groups,nbFields);
- K.scale(0.);
- K.axpy(sol);
- Unp.scale(0.);
- Unp.axpy(sol);
-
- for(int j=0; j<orderRK;j++){
- if(j){
- K.scale(b[j]);
- K.axpy(sol);
- }
-
- if (limiter) limiter->apply(K, groups);
- this->residual(claw,groups,K,resd);
- K.scale(0.);
- for(int k=0; k < groups.getNbElementGroups(); k++) {
- dgGroupOfElements *group = groups.getElementGroup(k);
- int nbNodes = group->getNbNodes();
- for(int i=0;i<group->getNbElements();i++) {
- residuEl.setAsProxy(resd.getGroupProxy(k),i*nbFields,nbFields);
- KEl.setAsProxy(K.getGroupProxy(k),i*nbFields,nbFields);
- iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
- iMassEl.mult(residuEl,KEl);
- }
- }
- Unp.axpy(K,a[j]);
- }
- if (limiter) limiter->apply(Unp, groups);
- sol.scale(0.);
- sol.axpy(Unp);
-}
-
-void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
+/*void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conservation law
dgGroupCollection &groups,
double h, // time-step
dgDofContainer &sol,
@@ -397,16 +339,16 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
fullMatrix<double> residuEl, KEl;
fullMatrix<double> iMassEl;
- int nbFields = claw.nbFields();
+ int nbFields = claw.getNbFields();
dgDofContainer **K;
K=new dgDofContainer*[nStages];
for(int i=0;i<nStages;i++){
- K[i]=new dgDofContainer(groups,nbFields);
+ K[i]=new dgDofContainer(&groups,nbFields);
K[i]->scale(0.);
}
- dgDofContainer Unp (groups,nbFields);
- dgDofContainer tmp (groups,nbFields);
+ dgDofContainer Unp (&groups,nbFields);
+ dgDofContainer tmp (&groups,nbFields);
Unp.scale(0.0);
Unp.axpy(sol);
@@ -435,7 +377,7 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
}
}
}
- this->residual(claw,groups,tmp,resd);
+ dgAlgorithm::residual(claw,groups,tmp,resd);
for(int k=0;k < groups.getNbElementGroups(); k++) {
if (k==0) {
b=b2;
@@ -464,7 +406,7 @@ void dgAlgorithm::multirateRungeKutta (const dgConservationLaw &claw, // conse
delete K[i];
}
delete []K;
- }
+ }*/
void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
@@ -474,7 +416,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
fullMatrix<double> &residual // residual !! at faces nodes
)
{
- int nbFields = claw.nbFields();
+ int nbFields = claw.getNbFields();
int nbNodesLeft = group.getGroupLeft().getNbNodes();
const dgBoundaryCondition *boundaryCondition = claw.getBoundaryCondition(group.getBoundaryTag());
// ----- 1 ---- get the solution at quadrature points
@@ -559,7 +501,7 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
dgDofContainer &residu)
{
solution.scatter();
- int nbFields=claw.nbFields();
+ int nbFields=claw.getNbFields();
//volume term
for(size_t i=0;i<groups.getNbElementGroups() ; i++) {
residu.getGroupProxy(i).scale(0);
@@ -618,7 +560,7 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
dataCacheDouble *maxConvectiveSpeed = claw.newMaxConvectiveSpeed(cacheMap);
dataCacheDouble *maximumDiffusivity = claw.newMaximumDiffusivity(cacheMap);
- const int nbFields = claw.nbFields();
+ const int nbFields = claw.getNbFields();
/* This is an estimate on how lengths changes with p
It is merely the smallest distance between gauss
points at order p + 1 */
@@ -647,3 +589,4 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
DT[iElement] = 1./spectralRadius;
}
}
+
diff --git a/Solver/dgAlgorithm.h b/Solver/dgAlgorithm.h
index 2baaa65305..f970010eb8 100644
--- a/Solver/dgAlgorithm.h
+++ b/Solver/dgAlgorithm.h
@@ -15,6 +15,7 @@ class dgSystemOfEquations;
class dgAlgorithm {
public :
+ static void registerBindings(binding *b);
static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
@@ -34,43 +35,19 @@ class dgAlgorithm {
fullMatrix<double> &solutionLeft,
fullMatrix<double> &residual // residual !! at faces nodes
);
- // works only if there is only 1 group of element
static void residual( const dgConservationLaw &claw, dgGroupCollection &groups,
dgDofContainer &solution, dgDofContainer &residual);
- void rungeKutta (
- const dgConservationLaw &claw,
- dgGroupCollection &groups,
- double h,
- dgDofContainer &residu,
- dgDofContainer &sol,
- dgLimiter *limiter=NULL,
- int orderRK=4);
static void computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group, // the group
fullMatrix<double> &solution, // the solution
std::vector<double> & DT // elementary time steps
);
- // works only (for the moment) for two groups of elements, small and big step
- // Uses RK43 from Schlegel (10 stages)
- void multirateRungeKutta (
- const dgConservationLaw &claw,
- dgGroupCollection &groups,
- double h,
- dgDofContainer &residu,
- dgDofContainer &sol,
- int orderRK=3);
static void multAddInverseMassMatrix ( /*dofManager &dof,*/
const dgGroupOfElements & group,
fullMatrix<double> &residu,
fullMatrix<double> &sol);
- static void buildGroups(GModel *model, int dim, int order,
- std::vector<dgGroupOfElements*> &eGroups,
- std::vector<dgGroupOfFaces*> &fGroups,
- std::vector<dgGroupOfFaces*> &bGroups,
- std::vector<dgGroupOfElements*> &ghostGroups);
};
-
#endif
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index a16e236ebf..aee3df289f 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -10,7 +10,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dgConservationLaw *_claw;
public:
term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
+ dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()),
cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
solutionRight(cacheMapRight.provideData("Solution")),
outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
@@ -34,7 +34,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dataCacheDouble &outsideValue;
public:
dirichlet(dgConservationLaw *claw, dataCacheMap &cacheMap,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMap, cacheMap.getNbEvaluationPoints(),claw->nbFields()),
+ dataCacheDouble(cacheMap, cacheMap.getNbEvaluationPoints(),claw->getNbFields()),
outsideValue(cacheMap.get(outsideValueFunctionName,this)){}
void _eval () {
for(int i=0;i<_value.size1();i++)
@@ -50,7 +50,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
dgConservationLaw *_claw;
public:
maximumDiffusivity(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
+ dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()),
cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
solutionRight(cacheMapRight.provideData("Solution")),
outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
@@ -89,7 +89,7 @@ class dgBoundaryConditionNeumann : public dgBoundaryCondition {
dataCacheDouble &flux;
public:
term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string fluxFunctionName):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
+ dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()),
flux(cacheMapLeft.get(fluxFunctionName,this))
{}
void _eval() {
@@ -113,7 +113,7 @@ class dgBoundarySymmetry : public dgBoundaryCondition {
dgConservationLaw *_claw;
public:
term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft):
- dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()), _claw(claw)
+ dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()), _claw(claw)
{
riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapLeft);
riemannSolver->addMeAsDependencyOf(this);
@@ -138,7 +138,7 @@ class dgBoundaryCondition0Flux : public dgBoundaryCondition {
class term : public dataCacheDouble {
public:
term(dgConservationLaw *claw,dataCacheMap &cacheMapLeft):
- dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()) {}
+ dataCacheDouble(cacheMapLeft,cacheMapLeft.getNbEvaluationPoints(),claw->getNbFields()) {}
void _eval() {
}
};
@@ -238,6 +238,8 @@ void dgConservationLaw::registerBindings(binding *b){
cm = cb->addMethod("newNeumannBoundary",&dgConservationLaw::newNeumannBoundary);
cm->setDescription("Create a new boundary condition with a given flux (no other fluxes will be computed, nor with the rieman solver nor the IP diffusive term");
cm->setArgNames("flux",NULL);
+ cm = cb->addMethod("getNbFields",&dgConservationLaw::getNbFields);
+ cm->setDescription("Return the number of fields composing the unknowns of this conservation law. For vectorial fields, each components is counted as one field.");
}
void dgBoundaryCondition::registerBindings(binding *b){
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 5013dff493..8058c567ac 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -35,7 +35,7 @@ class dgConservationLaw {
public:
virtual ~dgConservationLaw () {}
- int nbFields() const {return _nbf;}
+ int getNbFields() const {return _nbf;}
virtual dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newDiffusiveFlux (dataCacheMap &cacheMap) const {return NULL;}
diff --git a/Solver/dgDofContainer.cpp b/Solver/dgDofContainer.cpp
index fe13a95aad..bb0056c342 100644
--- a/Solver/dgDofContainer.cpp
+++ b/Solver/dgDofContainer.cpp
@@ -1,13 +1,18 @@
#include "GmshConfig.h"
#include "dgDofContainer.h"
+#include "function.h"
#include "dgGroupOfElements.h"
+#include "dgConservationLaw.h"
+#include "dgAlgorithm.h"
#ifdef HAVE_MPI
#include "mpi.h"
#else
#include "string.h"
#endif
-dgDofContainer::dgDofContainer (dgGroupCollection &groups, int nbFields):
- _groups(groups)
+#include <sstream>
+#include "MElement.h"
+dgDofContainer::dgDofContainer (dgGroupCollection *groups, int nbFields):
+ _groups(*groups)
{
int _dataSize = 0;
_dataSizeGhost=0;
@@ -173,8 +178,8 @@ double dgDofContainer::norm() {
#endif
}
void dgDofContainer::save(const std::string name) {
- FILE *f = fopen (name.c_str(),"rb");
- _data->binaryLoad(f);
+ FILE *f = fopen (name.c_str(),"wb");
+ _data->binarySave(f);
fclose(f);
}
void dgDofContainer::load(const std::string name) {
@@ -182,3 +187,112 @@ void dgDofContainer::load(const std::string name) {
_data->binaryLoad(f);
fclose(f);
}
+
+void dgDofContainer::L2Projection(std::string functionName){
+ dgDofContainer rhs(&_groups, _nbFields);
+ for (int iGroup=0;iGroup<_groups.getNbElementGroups();iGroup++) {
+ const dgGroupOfElements &group = *_groups.getElementGroup(iGroup);
+ fullMatrix<double> Source = fullMatrix<double> (group.getNbIntegrationPoints(),group.getNbElements() * _nbFields);
+ dataCacheMap cacheMap(group.getNbIntegrationPoints());
+ dataCacheElement &cacheElement = cacheMap.getElement();
+ cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
+ dataCacheDouble &sourceTerm = cacheMap.get(functionName);
+ fullMatrix<double> source;
+ for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
+ cacheElement.set(group.getElement(iElement));
+ source.setAsProxy(Source, iElement*_nbFields, _nbFields);
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iElement, iPt);
+ for (int k=0;k<_nbFields;k++)
+ source(iPt,k) = sourceTerm(iPt,k)*detJ;
+ }
+ }
+ rhs.getGroupProxy(iGroup).gemm(group.getSourceRedistributionMatrix(), Source);
+ dgAlgorithm::multAddInverseMassMatrix(group, rhs.getGroupProxy(iGroup), getGroupProxy(iGroup));
+ }
+}
+
+
+void dgDofContainer::exportMsh(const std::string name){
+ // the elementnodedata::export does not work !!
+
+ for (int ICOMP = 0; ICOMP<_nbFields;++ICOMP){
+ std::ostringstream name_oss;
+ name_oss<<name<<"_COMP_"<<ICOMP<<".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",name.c_str());
+ 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);
+ 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 ",sol(k,ICOMP));
+ }
+ 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;
+ */
+}
+
+
+#include "LuaBindings.h"
+void dgDofContainer::registerBindings(binding *b){
+ classBinding *cb = b->addClass<dgDofContainer>("dgDofContainer");
+ cb->setDescription("The DofContainer class provides a vector containing the degree of freedoms");
+ methodBinding *cm;
+ cm = cb->setConstructor<dgDofContainer,dgGroupCollection*,int>();
+ cm->setDescription("Build a vector");
+ cm->setArgNames("GroupCollection","nbFields",NULL);
+ cm = cb->addMethod("L2Projection",&dgDofContainer::L2Projection);
+ cm->setDescription("Project a function onto this vector");
+ cm->setArgNames("functionName",NULL);
+ cm = cb->addMethod("exportMsh",&dgDofContainer::exportMsh);
+ cm->setDescription("Export the dof for gmsh visualization");
+ cm->setArgNames("filename",NULL);
+}
diff --git a/Solver/dgDofContainer.h b/Solver/dgDofContainer.h
index 6c21641777..bcbdb732ea 100644
--- a/Solver/dgDofContainer.h
+++ b/Solver/dgDofContainer.h
@@ -29,12 +29,17 @@ public:
void axpy(dgDofContainer &x, double a=1.);
inline fullMatrix<double> &getGroupProxy(int gId){ return *(_dataProxys[gId]); }
inline fullMatrix<double> &getGroupProxy(const dgGroupOfElements* g){ return *(_dataProxys[_groupId[g]]); }
- dgDofContainer (dgGroupCollection &groups, int nbFields);
+ dgDofContainer (dgGroupCollection *groups, int nbFields);
~dgDofContainer ();
int getNbElements() {return _totalNbElements;}
void scatter();
void save(const std::string name);
void load(const std::string name);
void setAll(double v);
+ void L2Projection(std::string functionName);
+ void exportMsh(const std::string name);
+
+ static void registerBindings(binding *b);
+ inline dgGroupCollection *getGroups() { return &_groups; }
};
#endif
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 96ada631f2..c1cfdab539 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -1,5 +1,9 @@
+#include <iostream>
#include "GmshConfig.h"
#include "dgGroupOfElements.h"
+#include "dgConservationLaw.h"
+#include "dgDofContainer.h"
+#include "dgAlgorithm.h"
#include "MElement.h"
#include "polynomialBasis.h"
#include "Numeric.h"
@@ -635,6 +639,7 @@ void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
//static void partitionGroup (std::vector<MElement *>,){
//}
+
// this function creates groups of elements
// First, groups of faces are created on every physical group
// of dimension equal to the dimension of the problem minus one
@@ -642,28 +647,66 @@ void dgAlgorithm::partitionGroup(dgGroupOfElements &eGroup,
// -) Elements of different types are separated
// -) Then those groups are split into partitions
// -) Finally, those groups are re-numbered
-// Finally, group of interfaces are created
-// -) Groups of faces internal to a given group
-// -) Groups of faces between groups.
-void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
+void dgGroupCollection::buildGroupsOfElements(GModel *model, int dim, int order)
{
+ if(_groupsOfElementsBuilt)
+ return;
+ _groupsOfElementsBuilt=true;
_model = model;
- std::map<const std::string,std::set<MVertex*> > boundaryVertices;
- std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
- std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
std::vector<GEntity*> entities;
model->getEntities(entities);
std::map<int, std::vector<MElement *> >localElements;
- std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
int nlocal=0;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity *entity = entities[i];
+ if(entity->dim() == dim){
+ for (int iel=0; iel<entity->getNumMeshElements(); iel++){
+ MElement *el=entity->getMeshElement(iel);
+ if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
+ localElements[el->getType()].push_back(el);
+ nlocal++;
+ }
+ }
+ }
+ }
+
+
+ Msg::Info("%d groups of elements",localElements.size());
+ // do a group of elements for every element type in the mesh
+ for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
+ _elementGroups.push_back(new dgGroupOfElements(it->second,order));
+ }
+}
+
+// Finally, group of interfaces are created
+// -) Groups of faces internal to a given group
+// -) Groups of faces between groups.
+void dgGroupCollection::buildGroupsOfInterfaces() {
+ if(_groupsOfInterfacesBuilt)
+ return;
+ _groupsOfInterfacesBuilt=true;
+
+ int dim = _elementGroups[0]->getElement(0)->getDim();
+ int order = _elementGroups[0]->getOrder();
+
+ std::map<const std::string,std::set<MVertex*> > boundaryVertices;
+ std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
+ std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
+
+ std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
int nghosts=0;
- std::multimap<MElement*, short> &ghostsCells = model->getGhostCells();
+ std::multimap<MElement*, short> &ghostsCells = _model->getGhostCells();
+
+ std::vector<GEntity*> entities;
+ _model->getEntities(entities);
+
+
for(unsigned int i = 0; i < entities.size(); i++){
GEntity *entity = entities[i];
if(entity->dim() == dim-1){
for(unsigned int j = 0; j < entity->physicals.size(); j++){
- const std::string physicalName = model->getPhysicalName(entity->dim(), entity->physicals[j]);
+ const std::string physicalName = _model->getPhysicalName(entity->dim(), entity->physicals[j]);
for (int k = 0; k < entity->getNumMeshElements(); k++) {
MElement *element = entity->getMeshElement(k);
switch(dim) {
@@ -681,13 +724,11 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
}
}
}
- }else if(entity->dim() == dim){
+ }
+ else if(entity->dim() == dim){
for (int iel=0; iel<entity->getNumMeshElements(); iel++){
MElement *el=entity->getMeshElement(iel);
- if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
- localElements[el->getType()].push_back(el);
- nlocal++;
- }else{
+ if( ! (el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0) ){
std::multimap<MElement*, short>::iterator ghost=ghostsCells.lower_bound(el);
while(ghost!= ghostsCells.end() && ghost->first==el){
if(abs(ghost->second)-1==Msg::GetCommRank()){
@@ -700,13 +741,6 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
}
}
}
-
-
- Msg::Info("%d groups of elements",localElements.size());
- // do a group of elements for every element type in the mesh
- for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
- _elementGroups.push_back(new dgGroupOfElements(it->second,order));
- }
for (int i=0;i<_elementGroups.size();i++){
@@ -839,244 +873,40 @@ void dgGroupCollection::buildGroups(GModel *model, int dim, int order)
}
-void dgGroupCollection::buildGroups(GModel *model, int dim, int order,std::string groupType)
-{
- _model = model;
- std::map<const std::string,std::set<MVertex*> > boundaryVertices;
- std::map<const std::string,std::set<MEdge, Less_Edge> > boundaryEdges;
- std::map<const std::string,std::set<MFace, Less_Face> > boundaryFaces;
- std::vector<GEntity*> entities;
- model->getEntities(entities);
- std::map<int, std::vector<MElement *> >localElements;
- std::vector<std::map<int, std::vector<MElement *> > >ghostElements(Msg::GetCommSize()); // [partition][elementType]
- int nlocal=0;
- int nghosts=0;
- std::multimap<MElement*, short> &ghostsCells = model->getGhostCells();
- for(unsigned int i = 0; i < entities.size(); i++){
- GEntity *entity = entities[i];
- if(entity->dim() == dim-1){
- for(unsigned int j = 0; j < entity->physicals.size(); j++){
- const std::string physicalName = model->getPhysicalName(entity->dim(), entity->physicals[j]);
- for (int k = 0; k < entity->getNumMeshElements(); k++) {
- MElement *element = entity->getMeshElement(k);
- switch(dim) {
- case 1:
- boundaryVertices[physicalName].insert( element->getVertex(0) );
- break;
- case 2:
- boundaryEdges[physicalName].insert( element->getEdge(0) );
- break;
- case 3:
- boundaryFaces[physicalName].insert( element->getFace(0));
- break;
- default :
- throw;
- }
- }
- }
- }else if(entity->dim() == dim){
- double max_edgeM=-1.e22;
- double max_edgem=1.e22;
- double min_edgeM=-1.e22;
- double min_edgem=1.e22;
-
- for (int iel=0; iel<entity->getNumMeshElements(); iel++){
- MElement *el=entity->getMeshElement(iel);
- if (el->maxEdge()>max_edgeM) {
- max_edgeM=el->maxEdge();
- }
- if (el->maxEdge()<max_edgem) {
- max_edgem=el->maxEdge();
- }
- if (el->minEdge()>min_edgeM) {
- min_edgeM=el->minEdge();
- }
- if (el->minEdge()<min_edgeM) {
- min_edgem=el->minEdge();
- }
- }
- //printf("\n a1=%f, a2=%f \n",(min_edgeM+min_edgem)/2.,(max_edgeM+max_edgem)/2.);
- for (int iel=0; iel<entity->getNumMeshElements(); iel++){
- MElement *el=entity->getMeshElement(iel);
- if(el->getPartition()==Msg::GetCommRank()+1 || el->getPartition()==0){
-
- if (groupType=="minEdge") {
- if (el->minEdge()> (min_edgeM+min_edgem)/2.) {
- localElements[0].push_back(el);
- }
- else{
- localElements[1].push_back(el);
- }
- }
- else if(groupType=="maxEdge"){
- if (el->maxEdge()> (max_edgeM+max_edgem)/2.) {
- localElements[0].push_back(el);
- }
- else{
- localElements[1].push_back(el);
- }
-
- }
- else{
- localElements[el->getType()].push_back(el);
- }
+// Split the groups of elements depending on their local time step
+void dgGroupCollection::splitGroupsForMultirate(double refDt,dgConservationLaw *claw){
+ // 1) find the element with the smallest stable time step
+ double sr = DBL_MAX;
+ dgDofContainer *solution = new dgDofContainer((this),claw->getNbFields());
+ solution->scale(0.);
+ for (int i=0;i<getNbElementGroups();i++){
+ std::vector<double> DTS;
+ dgAlgorithm::computeElementaryTimeSteps(*claw, *getElementGroup(i), solution->getGroupProxy(i), DTS);
+ for (int k=0;k<DTS.size();k++){
+ std::cout << "SR: " << DTS[k] << std::endl;
+ sr = std::min(sr,DTS[k]);
+ }
+ }
+ delete solution;
+ #ifdef HAVE_MPI
+ double sr_min;
+ MPI_Allreduce((void *)&sr, &sr_min, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
+ sr=sr_min;
+ #endif
+ // sr is the time step for the most constrained element.
+ std::cout << "SR: " << sr << std::endl;
-
- nlocal++;
- }else{
- std::multimap<MElement*, short>::iterator ghost=ghostsCells.lower_bound(el);
- while(ghost!= ghostsCells.end() && ghost->first==el){
- if(abs(ghost->second)-1==Msg::GetCommRank()){
- ghostElements[el->getPartition()-1][el->getType()].push_back(el);
- nghosts+=1;
- }
- ghost++;
- }
- }
- }
- }
- }
-
-
- Msg::Info("%d groups of elements",localElements.size());
- for (int n=0; n<localElements.size(); n++) {
- printf("\n **** Group %i: %lu elements **** \n", n,localElements[n].size());
- }
- // do a group of elements for every element type in the mesh
- for (std::map<int, std::vector<MElement *> >::iterator it = localElements.begin(); it !=localElements.end() ; ++it){
- _elementGroups.push_back(new dgGroupOfElements(it->second,order));
- }
-
-
- for (int i=0;i<_elementGroups.size();i++){
- // create a group of faces on the boundary for every group ef elements
- switch(dim) {
- case 1 : {
- std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
- for(mapIt=boundaryVertices.begin(); mapIt!=boundaryVertices.end(); mapIt++) {
- dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
- if (gof->getNbElements())
- _boundaryGroups.push_back(gof);
- else
- delete gof;
- break;
- }
- }
- case 2 : {
- std::map<const std::string, std::set<MEdge, Less_Edge> >::iterator mapIt;
- for(mapIt=boundaryEdges.begin(); mapIt!=boundaryEdges.end(); mapIt++) {
- dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
- if(gof->getNbElements())
- _boundaryGroups.push_back(gof);
- else
- delete gof;
- }
- break;
- }
- case 3 : {
- std::map<const std::string, std::set<MFace, Less_Face> >::iterator mapIt;
- for(mapIt=boundaryFaces.begin(); mapIt!=boundaryFaces.end(); mapIt++) {
- dgGroupOfFaces *gof=new dgGroupOfFaces(*_elementGroups[i],mapIt->first,order,mapIt->second);
- if(gof->getNbElements())
- _boundaryGroups.push_back(gof);
- else
- delete gof;
- }
- break;
- }
- }
- // create a group of faces for every face that is common to elements of the same group
- _faceGroups.push_back(new dgGroupOfFaces(*_elementGroups[i],order));
- // create a groupe of d
- for (int j=i+1;j<_elementGroups.size();j++){
- dgGroupOfFaces *gof = new dgGroupOfFaces(*_elementGroups[i],*_elementGroups[j],order);
- if (gof->getNbElements())
- _faceGroups.push_back(gof);
- else
- delete gof;
- }
- }
- //create ghost groups
- for(int i=0;i<Msg::GetCommSize();i++){
- for (std::map<int, std::vector<MElement *> >::iterator it = ghostElements[i].begin(); it !=ghostElements[i].end() ; ++it){
- _ghostGroups.push_back(new dgGroupOfElements(it->second,order,i));
- }
- }
- //create face group for ghostGroups
- for (int i=0; i<_ghostGroups.size(); i++){
- for (int j=0;j<_elementGroups.size();j++){
- dgGroupOfFaces *gof = new dgGroupOfFaces(*_ghostGroups[i],*_elementGroups[j],order);
- if (gof->getNbElements())
- _faceGroups.push_back(gof);
- else
- delete gof;
- }
- }
-
- // build the ghosts structure
- int *nGhostElements = new int[Msg::GetCommSize()];
- int *nParentElements = new int[Msg::GetCommSize()];
- for (int i=0;i<Msg::GetCommSize();i++) {
- nGhostElements[i]=0;
- }
- for (size_t i = 0; i< getNbGhostGroups(); i++){
- dgGroupOfElements *group = getGhostGroup(i);
- nGhostElements[group->getGhostPartition()] += group->getNbElements();
- }
-#ifdef HAVE_MPI
- MPI_Alltoall(nGhostElements,1,MPI_INT,nParentElements,1,MPI_INT,MPI_COMM_WORLD);
-#else
- nParentElements[0]=nGhostElements[0];
-#endif
- int *shiftSend = new int[Msg::GetCommSize()];
- int *shiftRecv = new int[Msg::GetCommSize()];
- int *curShiftSend = new int[Msg::GetCommSize()];
- int totalSend=0,totalRecv=0;
- for(int i= 0; i<Msg::GetCommSize();i++){
- shiftSend[i] = (i==0 ? 0 : shiftSend[i-1]+nGhostElements[i-1]);
- curShiftSend[i] = shiftSend[i];
- shiftRecv[i] = (i==0 ? 0 : shiftRecv[i-1]+nParentElements[i-1]);
- totalSend += nGhostElements[i];
- totalRecv += nParentElements[i];
- }
-
- int *idSend = new int[totalSend];
- int *idRecv = new int[totalRecv];
- for (int i = 0; i<getNbGhostGroups(); i++){
- dgGroupOfElements *group = getGhostGroup(i);
- int part = group->getGhostPartition();
- for (int j=0; j< group->getNbElements() ; j++){
- idSend[curShiftSend[part]++] = group->getElement(j)->getNum();
- }
- }
-#ifdef HAVE_MPI
- MPI_Alltoallv(idSend,nGhostElements,shiftSend,MPI_INT,idRecv,nParentElements,shiftRecv,MPI_INT,MPI_COMM_WORLD);
-#else
- memcpy(idRecv,idSend,nParentElements[0]*sizeof(int));
-#endif
- //create a Map elementNum :: group, position in group
- std::map<int, std::pair<int,int> > elementMap;
- for(size_t i = 0; i< getNbElementGroups(); i++) {
- dgGroupOfElements *group = getElementGroup(i);
- for(int j=0; j< group->getNbElements(); j++){
- elementMap[group->getElement(j)->getNum()]=std::pair<int,int>(i,j);
- }
- }
- _elementsToSend.resize(Msg::GetCommSize());
- for(int i = 0; i<Msg::GetCommSize();i++){
- _elementsToSend[i].resize(nParentElements[i]);
- for(int j=0;j<nParentElements[i];j++){
- int num = idRecv[shiftRecv[i]+j];
- _elementsToSend[i][j]=elementMap[num];
- }
- }
- delete []idRecv;
- delete []idSend;
- delete []curShiftSend;
- delete []shiftSend;
- delete []shiftRecv;
}
+dgGroupCollection::dgGroupCollection()
+{
+ _groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
+}
+dgGroupCollection::dgGroupCollection(GModel *model, int dimension, int order)
+{
+ _groupsOfElementsBuilt=false;_groupsOfInterfacesBuilt=false;
+ buildGroupsOfElements(model,dimension,order);
+}
@@ -1090,3 +920,18 @@ dgGroupCollection::~dgGroupCollection() {
for (int i=0; i< _ghostGroups.size(); i++)
delete _ghostGroups[i];
}
+
+#include "LuaBindings.h"
+void dgGroupCollection::registerBindings(binding *b){
+ classBinding *cb = b->addClass<dgGroupCollection>("dgGroupCollection");
+ cb->setDescription("The GroupCollection class let you access to group partitioning functions");
+ methodBinding *cm;
+ cm = cb->setConstructor<dgGroupCollection,GModel*,int,int>();
+ cm->setDescription("Build the group of elements, separated by element type and element order. Additional partitioning can be done using splitGroupsForXXX specific functions");
+ cm->setArgNames("model","dimension","order",NULL);
+ cm = cb->addMethod("buildGroupsOfInterfaces",&dgGroupCollection::buildGroupsOfInterfaces);
+ 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("refDt","claw",NULL);
+}
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 3452aa31d4..063b0afd6e 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -20,6 +20,10 @@ class polynomialBasis;
class GEntity;
class GModel;
+class binding;
+class dgConservationLaw;
+
+
class dgElement {
MElement *_element;
// solution at points
@@ -192,7 +196,8 @@ class dgGroupCollection {
//{group,id} of the elements to send to each partition for a scatter operation
std::vector< std::vector<std::pair<int,int> > >_elementsToSend;
-
+ bool _groupsOfElementsBuilt;
+ bool _groupsOfInterfacesBuilt;
public:
inline int getNbElementGroups() const {return _elementGroups.size();}
inline int getNbFaceGroups() const {return _faceGroups.size();}
@@ -207,8 +212,14 @@ class dgGroupCollection {
inline int getImageElementGroup(int partId, int i) const {return _elementsToSend[partId][i].first;}
inline int getImageElementPositionInGroup(int partId, int i) const {return _elementsToSend[partId][i].second;}
- void buildGroups (GModel *model,int dimension, int order);
- void buildGroups (GModel *model,int dimension, int order, std::string groupType);
+ void buildGroupsOfElements (GModel *model,int dimension, int order);
+ void buildGroupsOfInterfaces ();
+
+ void splitGroupsForMultirate(double refDt,dgConservationLaw *claw);
+
+ dgGroupCollection(GModel *model, int dimension, int order);
+ dgGroupCollection();
+ static void registerBindings(binding *b);
~dgGroupCollection();
};
#endif
diff --git a/Solver/dgRungeKutta.cpp b/Solver/dgRungeKutta.cpp
new file mode 100644
index 0000000000..b07225b1a5
--- /dev/null
+++ b/Solver/dgRungeKutta.cpp
@@ -0,0 +1,90 @@
+#include "dgRungeKutta.h"
+#include "dgConservationLaw.h"
+#include "dgDofContainer.h"
+#include "dgLimiter.h"
+#include "dgAlgorithm.h"
+#include "dgGroupOfElements.h"
+#include "bindings.h"
+
+double dgRungeKutta::iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ double A[] = {0};
+ double b[] = {1};
+ return diagonalRK(claw, dt, solution, 1, A, b);
+}
+
+double dgRungeKutta::iterate44(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ double A[] = {0, 0.5, 0.5, 1};
+ double b[] = {1./6, 1./3, 1./3, 1./6};
+ return diagonalRK(claw, dt, solution, 4, A, b);
+}
+
+double dgRungeKutta::iterate22(const dgConservationLaw *claw, double dt, dgDofContainer *solution) {
+ double A[] = {0, 1};
+ double b[] = {1./2, 1./2};
+ return diagonalRK(claw, dt, solution, 2, A, b);
+}
+
+double dgRungeKutta::diagonalRK (const dgConservationLaw *claw,
+ double dt, // time-step
+ dgDofContainer *sol,
+ int nStages, double *A, double *b) {
+ // U_{n+1}=U_n+h*(SUM_i dt*b_i*K_i)
+ // K_i=M^(-1)R(U_n+dt*A_i*K_{i-1})
+ fullMatrix<double> residuEl, KEl;
+ fullMatrix<double> iMassEl;
+
+ int nbFields = claw->getNbFields();
+ dgGroupCollection *groups = sol->getGroups();
+
+ dgDofContainer K (groups, nbFields);
+ dgDofContainer Unp (groups, nbFields);
+ dgDofContainer resd(groups, nbFields);
+ K.scale(0.);
+ K.axpy(*sol);
+ Unp.scale(0.);
+ Unp.axpy(*sol);
+
+ for(int j=0; j<nStages;j++){
+ if(j){
+ K.scale(A[j]*dt);
+ K.axpy(*sol);
+ }
+
+ if (_limiter) _limiter->apply(K, *groups);
+ dgAlgorithm::residual(*claw,*groups,K,resd);
+ K.scale(0.);
+ for(int k=0; k < groups->getNbElementGroups(); k++) {
+ dgGroupOfElements *group = groups->getElementGroup(k);
+ int nbNodes = group->getNbNodes();
+ for(int i=0;i<group->getNbElements();i++) {
+ residuEl.setAsProxy(resd.getGroupProxy(k),i*nbFields,nbFields);
+ KEl.setAsProxy(K.getGroupProxy(k),i*nbFields,nbFields);
+ iMassEl.setAsProxy(group->getInverseMassMatrix(),i*nbNodes,nbNodes);
+ iMassEl.mult(residuEl,KEl);
+ }
+ }
+ Unp.axpy(K,b[j]*dt);
+ }
+ if (_limiter) _limiter->apply(Unp, *groups);
+ sol->scale(0.);
+ sol->axpy(Unp);
+}
+
+dgRungeKutta::dgRungeKutta():_limiter(NULL) {}
+
+void dgRungeKutta::registerBindings(binding *b) {
+ classBinding *cb = b->addClass<dgRungeKutta>("dgRungeKutta");
+ cb->setDescription("various explicit Runge-Kutta time stepping schemes");
+ methodBinding *cm;
+ cm = cb->setConstructor<dgRungeKutta>();
+ cm->setDescription("A new explicit runge kutta, pass parameters to the iterate function");
+ cm = cb->addMethod("iterateEuler",&dgRungeKutta::iterateEuler);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing a forward euler step of time step dt for the conservation law");
+ cm = cb->addMethod("iterate22",&dgRungeKutta::iterate22);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing Heun second order runge-kutta step of time step dt for the conservation law");
+ cm = cb->addMethod("iterate44",&dgRungeKutta::iterate44);
+ cm->setArgNames("law","dt","solution",NULL);
+ cm->setDescription("update solution by doing classical fourth order runge-kutta step of time step dt for the conservation law");
+}
diff --git a/Solver/dgRungeKutta.h b/Solver/dgRungeKutta.h
new file mode 100644
index 0000000000..c4a2495dc9
--- /dev/null
+++ b/Solver/dgRungeKutta.h
@@ -0,0 +1,19 @@
+#ifndef _DG_RUNGE_KUTTA_H_
+#define _DG_RUNGE_KUTTA_H_
+class dgConservationLaw;
+class dgDofContainer;
+class dgLimiter;
+class binding;
+class dgRungeKutta {
+ double diagonalRK(const dgConservationLaw *claw, double dt, dgDofContainer *solution, int nStages, double *A, double *b); // c == A
+ //static double nonDiagonalRK();
+ dgLimiter *_limiter;
+ public:
+ void setLimiter(dgLimiter *limiter) { _limiter = limiter; }
+ double iterateEuler(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
+ double iterate22(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
+ double iterate44(const dgConservationLaw *claw, double dt, dgDofContainer *solution);
+ static void registerBindings (binding *b);
+ dgRungeKutta();
+};
+#endif
diff --git a/Solver/dgSlopeLimiter.cpp b/Solver/dgSlopeLimiter.cpp
index f7c0ae70d1..cdef092477 100644
--- a/Solver/dgSlopeLimiter.cpp
+++ b/Solver/dgSlopeLimiter.cpp
@@ -6,14 +6,13 @@
//----------------------------------------------------------------------------------
bool dgSlopeLimiter::apply ( dgDofContainer &solution, dgGroupCollection &groups)
{
- printf("limit \n");
solution.scatter();
- int nbFields =_claw->nbFields();
+ int nbFields =_claw->getNbFields();
// first compute max and min of all fields for all stencils
//----------------------------------------------------------
- dgDofContainer MIN(groups, nbFields);
- dgDofContainer MAX(groups, nbFields);
+ dgDofContainer MIN(&groups, nbFields);
+ dgDofContainer MAX(&groups, nbFields);
MIN.setAll ( 1.e22);
MAX.setAll (-1.e22);
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index abfd8a151c..e9fa0c2e58 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -2,28 +2,16 @@
#include <stdlib.h>
#include <sstream>
#include "dgSystemOfEquations.h"
+#include "dgAlgorithm.h"
#include "function.h"
#include "MElement.h"
#include "PView.h"
#include "PViewData.h"
#include "dgLimiter.h"
+#include "dgRungeKutta.h"
#ifdef HAVE_MPI
#include "mpi.h"
#endif
-
-class dgConservationLawL2Projection : public dgConservationLaw {
- std::string _functionName;
-public:
- dgConservationLawL2Projection(const std::string & functionName, dgConservationLaw &_claw) :
- _functionName(functionName)
- {
- _nbf =_claw.nbFields();
- }
- dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
- return &cacheMap.get(_functionName);
- }
-};
-
dgSystemOfEquations::dgSystemOfEquations(GModel *gm){
_gm=gm;
_dimension = _gm->getNumRegions() ? 3 : _gm->getNumFaces() ? 2 : 1 ;
@@ -46,6 +34,7 @@ static dgSystemOfEquations *myConstructorPtr(GModel* gm){
return new dgSystemOfEquations(gm);
}
+
void dgSystemOfEquations::registerBindings(binding *b) {
classBinding *cb = b->addClass<dgSystemOfEquations>("dgSystemOfEquations");
cb->setDescription("a class to rule them all :-) -- bad description, this class will be removed anyway");
@@ -58,9 +47,6 @@ void dgSystemOfEquations::registerBindings(binding *b) {
cm->setDescription("set the conservation law this system solves");
cm = cb->addMethod("setup",&dgSystemOfEquations::setup);
cm->setDescription("allocate and init internal stuff, call this function after setOrder and setLaw and before anything else on this instance");
- cm = cb->addMethod("createGroups",&dgSystemOfEquations::createGroups);
- cm->setArgNames("groupType",NULL);
- cm->setDescription("allocate and init internal stuff, creates groups form criterion groupType");
cm = cb->addMethod("exportSolution",&dgSystemOfEquations::exportSolution);
cm->setArgNames("filename",NULL);
cm->setDescription("Print the solution into a file. This file does not contain the mesh. To visualize the solution in gmsh you have to open the mesh file first.");
@@ -83,9 +69,9 @@ void dgSystemOfEquations::registerBindings(binding *b) {
cm = cb->addMethod("RK44_limiter",&dgSystemOfEquations::RK44_limiter);
cm->setArgNames("dt",NULL);
cm->setDescription("do one RK44 time step with the slope limiter (only for p=1)");
- cm = cb->addMethod("multirateRK43",&dgSystemOfEquations::multirateRK43);
- cm->setArgNames("dt",NULL);
- cm->setDescription("Do a runge-kuta temporal iteration with 4 sub-steps and a precision order of 3 using different time-step depending on the element size. This function returns the sum of the nodal residuals.");
+ ////cm = cb->addMethod("multirateRK43",&dgSystemOfEquations::multirateRK43);
+ //cm->setArgNames("dt",NULL);
+ //cm->setDescription("Do a runge-kuta temporal iteration with 4 sub-steps and a precision order of 3 using different time-step depending on the element size. This function returns the sum of the nodal residuals.");
cm = cb->addMethod("saveSolution",&dgSystemOfEquations::saveSolution);
cm->setArgNames("filename",NULL);
cm->setDescription("dump the solution in binary format");
@@ -96,35 +82,22 @@ void dgSystemOfEquations::registerBindings(binding *b) {
// do a L2 projection
void dgSystemOfEquations::L2Projection (std::string functionName){
- dgConservationLawL2Projection Law(functionName,*_claw);
- for (int i=0;i<_groups.getNbElementGroups();i++){
- dgGroupOfElements *group = _groups.getElementGroup(i);
- _algo->residualVolume(Law,*group,_solution->getGroupProxy(i),_rightHandSide->getGroupProxy(i));
- _algo->multAddInverseMassMatrix(*group,_rightHandSide->getGroupProxy(i),_solution->getGroupProxy(i));
- }
+ _solution->L2Projection(functionName);
}
-// dgSystemOfEquations::setup() + build groups with criterion:
-// - default: elementType
-// - minEdge (based on minimum edges of elements) , for the moment only two groups possible
-// - maxEdge (based on maximum edges of elements) , for the moment only two groups possible
-void dgSystemOfEquations::createGroups(std::string groupType){
- _groups.buildGroups(_gm,_dimension,_order,groupType);
- _solution = new dgDofContainer(_groups,_claw->nbFields());
- _rightHandSide = new dgDofContainer(_groups,_claw->nbFields());
- }
-
// ok, we can setup the groups and create solution vectors
void dgSystemOfEquations::setup(){
if (!_claw) throw;
- _groups.buildGroups(_gm,_dimension,_order);
- _solution = new dgDofContainer(_groups,_claw->nbFields());
- _rightHandSide = new dgDofContainer(_groups,_claw->nbFields());
+ _groups.buildGroupsOfElements(_gm,_dimension,_order);
+ _groups.buildGroupsOfInterfaces();
+ _solution = new dgDofContainer(&_groups,_claw->getNbFields());
+ _rightHandSide = new dgDofContainer(&_groups,_claw->getNbFields());
}
double dgSystemOfEquations::RK44(double dt){
- _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide);
+ dgRungeKutta rk;
+ rk.iterate44(_claw, dt, _solution);
return _solution->norm();
}
@@ -132,7 +105,7 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
double sr = 1.e22;
for (int i=0;i<_groups.getNbElementGroups();i++){
std::vector<double> DTS;
- _algo->computeElementaryTimeSteps(*_claw, *_groups.getElementGroup(i), _solution->getGroupProxy(i), DTS);
+ dgAlgorithm::computeElementaryTimeSteps(*_claw, *_groups.getElementGroup(i), _solution->getGroupProxy(i), DTS);
for (int k=0;k<DTS.size();k++) sr = std::min(sr,DTS[k]);
}
#ifdef HAVE_MPI
@@ -146,22 +119,25 @@ double dgSystemOfEquations::computeInvSpectralRadius(){
double dgSystemOfEquations::RK44_limiter(double dt){
dgLimiter *sl = new dgSlopeLimiter(_claw);
- _algo->rungeKutta(*_claw,_groups, dt, *_solution, *_rightHandSide, sl);
+ dgRungeKutta rk;
+ rk.setLimiter(sl);
+ rk.iterate44(_claw, dt, _solution);
delete sl;
return _solution->norm();
}
double dgSystemOfEquations::ForwardEuler(double dt){
- _algo->rungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide, NULL,1);
- return _solution->norm();
-}
-double dgSystemOfEquations::multirateRK43(double dt){
- _algo->multirateRungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide);
+ dgRungeKutta rk;
+ rk.iterateEuler(_claw, dt, _solution);
return _solution->norm();
}
+//double dgSystemOfEquations::multirateRK43(double dt){
+ //dgAlgorithm::multirateRungeKutta(*_claw, _groups, dt, *_solution, *_rightHandSide);
+ //return _solution->norm();
+//}
void dgSystemOfEquations::exportSolution(std::string outputFile){
- export_solution_as_is(outputFile);
+ _solution->exportMsh(outputFile);
}
void dgSystemOfEquations::limitSolution(){
@@ -185,72 +161,3 @@ void dgSystemOfEquations::saveSolution (std::string name) {
void dgSystemOfEquations::loadSolution (std::string name){
_solution->load(name);
}
-
-void dgSystemOfEquations::export_solution_as_is (const std::string &name){
- // the elementnodedata::export does not work !!
-
- for (int ICOMP = 0; ICOMP<_claw->nbFields();++ICOMP){
- std::ostringstream name_oss;
- name_oss<<name<<"_COMP_"<<ICOMP<<".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",name.c_str());
- 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);
- for (int iElement=0 ; iElement< group->getNbElements() ;++iElement) {
- MElement *e = group->getElement(iElement);
- int num = e->getNum();
- fullMatrix<double> sol = getSolutionProxy (i, iElement);
- fprintf(f,"%d %d",num,sol.size1());
- for (int k=0;k<sol.size1();++k) {
- fprintf(f," %.16E ",sol(k,ICOMP));
- }
- 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 = getSolutionProxy (i, iElement);
- 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;
-}
-
-
-
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index 85625e8eea..0cafe0ced7 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -4,9 +4,7 @@
#include <utility>
#include "GmshConfig.h"
#include "GModel.h"
-#include "dgAlgorithm.h"
#include "dgGroupOfElements.h"
-#include "dgAlgorithm.h"
#include "dgConservationLaw.h"
#include "Gmsh.h"
#include "dgLimiter.h"
@@ -21,8 +19,6 @@ class dgSystemOfEquations {
// the mesh and the model
GModel *_gm;
dgGroupCollection _groups;
- // the algorithm that computes DG operators
- dgAlgorithm *_algo;
// the conservation law
dgConservationLaw *_claw;
std::string _cLawName;
@@ -42,25 +38,21 @@ public:
void setOrder (int order); // set the polynomial order
void setConservationLaw (dgConservationLaw *law); // set the conservationLaw
dgSystemOfEquations(GModel *_gm);
- void createGroups(std::string groupType); // create groups from criterion: minEdge (2 groups), maxEdge (2 groups), elementType, allocate (same as setup())
void setup (); // setup the groups and allocate
void exportSolution (std::string filename); // export the solution
void limitSolution (); // apply the limiter on the solution
- double RK44 (double dt);
- double RK44_limiter (double dt);
- double ForwardEuler (double dt);
- double multirateRK43 (double dt);
void L2Projection (std::string functionName); // assign the solution to a given function
+ double RK44(double dt);
+ double RK44_limiter(double dt);
+ double ForwardEuler(double dt);
+
double computeInvSpectralRadius();
static void registerBindings(binding *b);
- inline const fullMatrix<double> getSolutionProxy (int iGroup, int iElement){
- return fullMatrix<double> ( _solution->getGroupProxy(iGroup), iElement * _claw->nbFields(),_claw->nbFields());
- }
- void export_solution_as_is (const std::string &fileName);
void saveSolution (std::string fileName) ;
void loadSolution (std::string fileName);
+
~dgSystemOfEquations();
};
diff --git a/contrib/kbipack/gmp_matrix_io.c b/contrib/kbipack/gmp_matrix_io.c
index f85639304a..ca1c6e4877 100644
--- a/contrib/kbipack/gmp_matrix_io.c
+++ b/contrib/kbipack/gmp_matrix_io.c
@@ -51,7 +51,7 @@ gmp_matrix * gmp_matrix_read_coord(char* filename)
}
/* First read the size of the matrix */
- read_values = sscanf(buffer, "%u %u %u", &nrows, &ncols, &dummy);
+ read_values = sscanf(buffer, "%lu %lu %lu", &nrows, &ncols, &dummy);
/* Create the matrix */
new_matrix = create_gmp_matrix_zero(nrows, ncols);
--
GitLab