From 127d0888600c9ef829ba3658d50d9baf1d4d0d5d Mon Sep 17 00:00:00 2001
From: Jonathan Lambrechts <jonathan.lambrechts@uclouvain.be>
Date: Tue, 2 Feb 2010 15:59:52 +0000
Subject: [PATCH] dg: add neumann bnd condition, add maxDiffusivityOutside on
dirichlet bnd
---
Solver/TESTCASES/supra.lua | 42 ++++++++++++++------
Solver/dgAlgorithm.cpp | 6 ++-
Solver/dgConservationLaw.cpp | 77 +++++++++++++++++++++++++++++++++++-
Solver/dgConservationLaw.h | 2 +
4 files changed, 113 insertions(+), 14 deletions(-)
diff --git a/Solver/TESTCASES/supra.lua b/Solver/TESTCASES/supra.lua
index 69539b7708..5b6a7e476c 100644
--- a/Solver/TESTCASES/supra.lua
+++ b/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
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index fde2d637db..4e4c426205 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/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);
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index 19eb442167..a16e236ebf 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/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){
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 196b9b8b51..5013dff493 100644
--- a/Solver/dgConservationLaw.h
+++ 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();
--
GitLab