Precomputed inner radius

Solver/TESTCASES/MultirateAdvection.lua

@@ -53,7 +53,7 @@ FS = functionLua(1, 'initial_condition', {'XYZ'}):getName()
 GC=dgGroupCollection(myModel,2,order)
 solTmp=dgDofContainer(GC,1)
 solTmp:L2Projection(FS)
-dt=GC:splitGroupsForMultirate(20,law,solTmp)
+dt=GC:splitGroupsForMultirate(100,law,solTmp)
 GC:buildGroupsOfInterfaces(myModel,2,order)
 solution=dgDofContainer(GC,1)
 solution:L2Projection(FS)
@@ -87,7 +87,7 @@ time=0
 -- multirateRK:setLimiter(limiter)
 --for i=1,1000
 i=0
-while time<0.1 do
+while time<0.2 do
 --     norm1 = RK:iterate43SchlegelJCAM2009(law,dt,solution)
 -- TEST with Explicit Euler multirate !!!
     norm2 = multirateRK:iterate(dt,solution2)

Solver/TESTCASES/rect.geo

@@ -18,7 +18,7 @@ Plane Surface(11) = {10};
 Physical Surface("sprut") = {11, 9};
 Physical Line("Walls") = {5, 7, 3, 4, 1};
 Physical Line("Top") = {6};
-//Recombine Surface {9, 11};
+Recombine Surface {9, 11};
 
 Field[1] = MathEval;
 Field[1].F = "0.01*1+0.01*100*(y*y)";

Solver/dgAlgorithm.cpp

@@ -107,7 +107,7 @@ void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF man
     sol.setAsProxy(solution, iElement*nbFields, nbFields);
     MElement *e = group.getElement(iElement);
     cacheElement.set(e);
-    const double L = e->minEdge();
+    const double L = group.getInnerRadius(iElement);
     double spectralRadius = 0.0;
     if (maximumDiffusivity){
       double nu = (*maximumDiffusivity)(0,0);

Solver/dgGroupOfElements.cpp

@@ -74,12 +74,14 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
   _imass = new fullMatrix<double> (nbPsi,nbPsi*e.size()); 
   _dPsiDx = new fullMatrix<double> ( _integration->size1()*3,nbPsi*e.size());
   _elementVolume = new fullMatrix<double> (e.size(),1);
+  _innerRadii = new fullMatrix<double> (e.size(),1);
   double g[256][3],f[256];
   
   for (int i=0;i<_elements.size();i++){
     MElement *e = _elements[i];
     element_to_index[e] = i;
     fullMatrix<double> imass(*_imass,nbPsi*i,nbPsi);
+    (*_innerRadii)(i,0)=e->getInnerRadius();
     for (int j=0;j< _integration->size1() ; j++ ){
       _fs.f((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), f);
       _fs.df((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), g);
@@ -132,7 +134,9 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e,
 }
 
 void dgGroupOfElements::copyPrivateDataFrom(const dgGroupOfElements *from){
-//  Msg::Error("dgGroupOfElements::copyPrivateDataFrom not yet implemented, ask Richard");
+  _multirateExponent=from->getMultirateExponent();
+  _multirateInnerBuffer=from->getIsInnerMultirateBuffer();
+  _multirateOuterBuffer=from->getIsOuterMultirateBuffer();
 }
 
 dgGroupOfElements::~dgGroupOfElements(){
@@ -1051,7 +1055,9 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
   Msg::Info("Splitting Groups for multirate time stepping");
   maxLevels--;// Number becomes maximum id
 
+  // Interfaces are not built yet, so we use a "mini" structure to deduce neighboring information
   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;
 
   std::vector<std::vector<int> >newGroupIds;
@@ -1113,7 +1119,6 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
   _dtMaxExponent=dtMaxExponent;
   std::vector<MElement *>currentNewGroup;
   std::vector< dgGroupOfElements* >newGroups;// indexed by newGroupId
-  std::map<int,int>oldToNewGroupsMap;// oldGroupId to newGroupId
 
   int lowerLevelGroupIdStart=-1;
   int lowerLevelGroupIdEnd=-1;
@@ -1211,7 +1216,6 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
             newGroup->_multirateInnerBuffer=false;
             newGroups.resize(currentNewGroupId+1);
             newGroups[currentNewGroupId]=newGroup;
-            oldToNewGroupsMap[it->first]=currentNewGroupId;
             currentNewGroupId++;
           }
 
@@ -1227,7 +1231,6 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
             newGroup->_multirateInnerBuffer=true;
             newGroups.resize(currentNewGroupId+1);
             newGroups[currentNewGroupId]=newGroup;
-            oldToNewGroupsMap[it->first]=currentNewGroupId;
             currentNewGroupId++;
           }
         }
@@ -1255,7 +1258,6 @@ double dgGroupCollection::splitGroupsForMultirate(int maxLevels,dgConservationLa
           newGroup->_multirateInnerBuffer=false;
           newGroups.resize(currentNewGroupId+1);
           newGroups[currentNewGroupId]=newGroup;
-          oldToNewGroupsMap[it->first]=currentNewGroupId;
           currentNewGroupId++;
         }
         else

Solver/dgGroupOfElements.h

@@ -80,6 +80,8 @@ class dgGroupOfElements {
   //
   // volume/surface/length of the element (sum_qp w_i detJ_i)
   fullMatrix<double> *_elementVolume;
+  // inradius of the element, i.e. a lengthscale to compute the timestep
+  fullMatrix<double> *_innerRadii;
   // forbid the copy 
   //  dgGroupOfElements (const dgGroupOfElements &e, int order) {}
   //  dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
@@ -107,6 +109,7 @@ public:
   inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
   inline const fullMatrix<double> & getSourceRedistributionMatrix () const {return *_redistributionSource;}
   inline double getElementVolume (int iElement)const {return (*_elementVolume)(iElement,0);}
+  inline double getInnerRadius(int iElement)const {return (*_innerRadii)(iElement,0);}
   inline double getDetJ (int iElement, int iGaussPoint) const {return (*_mapping)(iElement, 10*iGaussPoint + 9);}
   inline double getInvJ (int iElement, int iGaussPoint, int i, int j) const {return (*_mapping)(iElement, 10*iGaussPoint + i + 3*j);}
   inline fullMatrix<double> & getDPsiDx() const { return *_dPsiDx;}