dg: add neumann bnd condition, add maxDiffusivityOutside on dirichlet bnd

Solver/TESTCASES/supra.lua

@@ -4,12 +4,27 @@ model:load ('square_supra.msh')
 dg = dgSystemOfEquations (model)
 dg:setOrder(1)
 
+c=1e-3
+n=20
+t=0
+dt=1e-8
+T=3e-4
 
 function diffusivity( sol , f )
-  n=20
-  c=1e-3
-  for i=0,sol:size1()-1 do
-    f:set (i, 0, sol:get(i,0)^(n-1)*n/c )
+  for i=0,f:size1()-1 do
+    f:set (i, 0, math.abs(sol:get(i,0))^(n-1)*n/c )
+  end
+end
+
+function valueRight(f)
+  for i=0,f:size1()-1 do
+    f:set(i,0,math.sin(2*math.pi*t/T))
+  end
+end
+
+function valueLeft(f)
+  for i=0,f:size1()-1 do
+    f:set(i,0,-math.sin(2*math.pi*t/T))
   end
 end
 
@@ -17,16 +32,20 @@ law = dgConservationLawAdvectionDiffusion('',functionLua(1,'diffusivity',{'Solut
 dg:setConservationLaw(law)
 
 -- boundary condition
+--[[
 outsideLeft=fullMatrix(1,1)
 outsideLeft:set(0,0,-1)
 outsideRight=fullMatrix(1,1)
 outsideRight:set(0,0,1)
+--]]
 law:addBoundaryCondition('Top',law:new0FluxBoundary())
 law:addBoundaryCondition('Bottom',law:new0FluxBoundary())
-leftFunction=functionConstant(outsideLeft):getName()
-rightFunction=functionConstant(outsideRight):getName()
-print(leftFunction)
-print(rightFunction)
+leftFunction=functionLua(1,'valueLeft',{}):getName()
+rightFunction=functionLua(1,'valueRight',{}):getName()
+--[[
+law:addBoundaryCondition('Left',law:newNeumannBoundary(leftFunction))
+law:addBoundaryCondition('Right',law:newNeumannBoundary(rightFunction))
+--]]
 law:addBoundaryCondition('Left',law:newOutsideValueBoundary(leftFunction))
 law:addBoundaryCondition('Right',law:newOutsideValueBoundary(rightFunction))
 
@@ -43,13 +62,14 @@ function initial_condition( xyz , f )
 end
 dg:L2Projection(functionLua(1,'initial_condition',{'XYZ'}):getName())
 
-dg:exportSolution('output/Diffusion_00000')
+dg:exportSolution('output/supra-00000')
 
 -- main loop
 for i=1,150000 do
-  norm = dg:RK44_limiter(0.000005)
+  norm = dg:RK44_limiter(dt)
+  t=t+dt
   if (i % 100 == 0) then 
     print('iter',i,norm)
-    dg:exportSolution(string.format("output/Diffusion-%05d", i)) 
+    dg:exportSolution(string.format("output/supra-%05d", i)) 
   end
 end

Solver/dgAlgorithm.cpp

@@ -500,6 +500,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
   // viscous
   dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
   dataCacheDouble *maximumDiffusivityLeft = claw.newMaximumDiffusivity(cacheMapLeft);
+  dataCacheDouble *maximumDiffusivityOut = boundaryCondition->newMaximumDiffusivity(cacheMapLeft);
   dataCacheDouble *neumann   = boundaryCondition->newDiffusiveNeumannBC(cacheMapLeft);
   dataCacheDouble *dirichlet = boundaryCondition->newDiffusiveDirichletBC(cacheMapLeft);
 
@@ -540,13 +541,14 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
         for (int k=0;k<nbFields;k++) { 
           double minl = group.getElementVolumeLeft(iFace)/group.getInterfaceSurface(iFace);
           double nu = (*maximumDiffusivityLeft)(iPt,0);
-	  double mu = (p+1)*(p+dim)/dim*nu/minl;
+          if(maximumDiffusivityOut)
+            nu = std::max(nu,(*maximumDiffusivityOut)(iPt,0));
+          double mu = (p+1)*(p+dim)/dim*nu/minl;
           double solutionJumpPenalty = (solutionQPLeft(iPt,k)-(*dirichlet)(iPt,k))/2*mu;
           normalFluxQP(iPt,k) -= ((*neumann)(iPt,k)+solutionJumpPenalty)*detJ;
         }
       }
     }    
-
   }
   // ----- 3 ---- do the redistribution at face nodes using BLAS3
   residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);

Solver/dgConservationLaw.cpp

@@ -17,7 +17,8 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
       _claw(claw)
     {
       riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapRight);
-      riemannSolver->addMeAsDependencyOf(this);
+      if(riemannSolver)
+        riemannSolver->addMeAsDependencyOf(this);
     }
 
     void _eval() {
@@ -29,6 +30,45 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
       }
     }
   };
+  class dirichlet : public dataCacheDouble {
+    dataCacheDouble &outsideValue;
+    public:
+    dirichlet(dgConservationLaw *claw, dataCacheMap &cacheMap,const std::string outsideValueFunctionName):
+      dataCacheDouble(cacheMap, cacheMap.getNbEvaluationPoints(),claw->nbFields()),
+      outsideValue(cacheMap.get(outsideValueFunctionName,this)){}
+    void _eval () { 
+      for(int i=0;i<_value.size1();i++)
+        for(int j=0;j<_value.size2();j++)
+          _value(i,j)=outsideValue(i,j);
+    }
+  };
+  class maximumDiffusivity : public dataCacheDouble {
+    dataCacheMap cacheMapRight; // new cacheMap to  pass to the Riemann solver
+    dataCacheDouble &solutionRight;
+    dataCacheDouble &outsideValue;
+    dataCacheDouble *maxDif;
+    dgConservationLaw *_claw;
+    public:
+    maximumDiffusivity(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
+      dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
+      cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
+      solutionRight(cacheMapRight.provideData("Solution")),
+      outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
+      _claw(claw)
+    {
+      maxDif = _claw->newMaximumDiffusivity(cacheMapRight);
+      if(maxDif)
+        maxDif->addMeAsDependencyOf(this);
+    }
+    void _eval() {
+      solutionRight.set(outsideValue());
+      if(maxDif){
+        for(int i=0;i<_value.size1(); i++)
+          for(int j=0;j<_value.size2(); j++)
+            _value(i,j) = (*maxDif)(i,j);
+      }
+    }
+  };
   public:
   dgBoundaryConditionOutsideValue(dgConservationLaw *claw,const std::string outsideValueFunctionName): dgBoundaryCondition(claw),
     _outsideValueFunctionName(outsideValueFunctionName)
@@ -36,6 +76,35 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
   dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
     return new term(_claw,cacheMapLeft,_outsideValueFunctionName);
   }
+  dataCacheDouble *newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const {
+    return new dirichlet(_claw,cacheMapLeft,_outsideValueFunctionName);
+  }
+  dataCacheDouble *newMaximumDiffusivity(dataCacheMap &cacheMapLeft) const {
+    return new maximumDiffusivity(_claw,cacheMapLeft,_outsideValueFunctionName); 
+  }
+};
+class dgBoundaryConditionNeumann : public dgBoundaryCondition {
+  std::string _fluxFunctionName;
+  class term : public dataCacheDouble {
+    dataCacheDouble &flux;
+    public:
+    term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string fluxFunctionName):
+      dataCacheDouble(cacheMapLeft, cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
+      flux(cacheMapLeft.get(fluxFunctionName,this))
+    {}
+    void _eval() {
+      for(int i=0;i<_value.size1(); i++)
+        for(int j=0;j<_value.size2(); j++)
+          _value(i,j) = (flux)(i,j);
+    }
+  };
+  public:
+  dgBoundaryConditionNeumann(dgConservationLaw *claw,const std::string fluxFunctionName): dgBoundaryCondition(claw),
+    _fluxFunctionName(fluxFunctionName)
+  { }
+  dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+    return new term(_claw,cacheMapLeft,_fluxFunctionName);
+  }
 };
 
 class dgBoundarySymmetry : public dgBoundaryCondition {
@@ -86,6 +155,9 @@ dgBoundaryCondition *dgConservationLaw::newSymmetryBoundary() {
 dgBoundaryCondition *dgConservationLaw::newOutsideValueBoundary(const std::string outsideValueFunctionName) {
   return new dgBoundaryConditionOutsideValue(this,outsideValueFunctionName);
 }
+dgBoundaryCondition *dgConservationLaw::newNeumannBoundary(const std::string fluxFunctionName) {
+  return new dgBoundaryConditionNeumann(this,fluxFunctionName);
+}
 dgBoundaryCondition *dgConservationLaw::new0FluxBoundary() {
   return new dgBoundaryCondition0Flux(this);
 }
@@ -163,6 +235,9 @@ void dgConservationLaw::registerBindings(binding *b){
   cm = cb->addMethod("newOutsideValueBoundary",&dgConservationLaw::newOutsideValueBoundary);
   cm->setDescription("Create a new boundary condition which compute the fluxes using the Riemann solver using the 'outsideFunction' function to compute external values.");
   cm->setArgNames("outsideFunction",NULL);
+  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);
 }
 
 void dgBoundaryCondition::registerBindings(binding *b){

Solver/dgConservationLaw.h

@@ -24,6 +24,7 @@ class dgBoundaryCondition {
   virtual dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const = 0;
   virtual dataCacheDouble *newDiffusiveNeumannBC(dataCacheMap &cacheMapLeft) const;
   virtual dataCacheDouble *newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const;
+  virtual dataCacheDouble *newMaximumDiffusivity(dataCacheMap &cacheMapLeft)const {return NULL;}
   static void registerBindings(binding *b);
 };
 
@@ -60,6 +61,7 @@ class dgConservationLaw {
 
   //a generic boundary condition using the Riemann solver of the conservation Law
   dgBoundaryCondition *newOutsideValueBoundary(std::string outsideValueFunctionName);
+  dgBoundaryCondition *newNeumannBoundary(const std::string fluxFunctionName);
   dgBoundaryCondition *new0FluxBoundary();
   dgBoundaryCondition *newSymmetryBoundary();