Skip to content
Snippets Groups Projects
Commit 4469c334 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

dg: IP

parent d25c8934
No related branches found
No related tags found
No related merge requests found
......@@ -289,6 +289,7 @@ bool fullMatrix<double>::svd(fullMatrix<double> &V, fullVector<double> &S)
return false;
}
#if defined(HAVE_LUA)
template<>
int fullMatrix<double>::gemm(lua_State *L)
......
......@@ -134,6 +134,7 @@ class fullMatrix
lua_pushnumber(L, (*this)(r,c));
return 1;
}
int set(lua_State *L){
int r = luaL_checkint(L, 1);
int c = luaL_checkint(L, 2);
......
......@@ -11,7 +11,11 @@ v=fullMatrix(3,1);
v:set(0,0,0.15)
v:set(1,0,0.05)
v:set(2,0,0)
dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v))
-- diffusivity
nu=fullMatrix(1,1);
nu:set(0,0,0.001)
dg:setConservationLaw('AdvectionDiffusion',createFunction.constant(v),createFunction.constant(nu))
-- boundary condition
outside=fullMatrix(1,1)
......
model = GModel ()
model:load ('square.geo')
model:load ('square.msh')
dg = dgSystemOfEquations (model)
dg:setOrder(5)
-- conservation law
-- advection speed
nu=fullMatrix(1,1);
nu:set(0,0,0.01)
dg:setConservationLaw('AdvectionDiffusion','',createFunction.constant(nu))
-- boundary condition
outside=fullMatrix(1,1)
outside:set(0,0,0.)
dg:addBoundaryCondition('Border','0Flux')
dg:setup()
-- initial condition
function initial_condition( _x , _f )
xyz = fullMatrix(_x)
f = fullMatrix(_f)
for i=0,xyz:size1()-1 do
x = xyz:get(i,0)
y = xyz:get(i,1)
z = xyz:get(i,2)
f:set (i, 0, math.exp(-100*((x-0.2)^2 +(y-0.3)^2)))
end
end
dg:L2Projection(createFunction.lua(1,'initial_condition','XYZ'))
dg:exportSolution('output/Diffusion_00000')
-- main loop
for i=1,10000 do
norm = dg:RK44(0.001)
if (i % 50 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Diffusion-%05d", i))
end
end
......@@ -16,7 +16,7 @@
void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
const fullMatrix<double> &solution,
fullMatrix<double> &solution,
fullMatrix<double> &residual // the residual
)
{
......@@ -45,27 +45,26 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
// provided dataCache
cacheMap.provideData("UVW").set(group.getIntegrationPointsMatrix());
dataCacheDouble &solutionQPe = cacheMap.provideData("Solution");
dataCacheDouble &gradSolutionQPe = cacheMap.provideData("gradSolution");
dataCacheDouble &gradientSolutionQPe = cacheMap.provideData("SolutionGradient");
gradientSolutionQPe.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
// needed dataCache
dataCacheDouble *sourceTerm = claw.newSourceTerm(cacheMap);
// fluxes in XYZ coordinates;
dataCacheDouble *convectiveFlux = claw.newConvectiveFlux(cacheMap);
dataCacheDouble *diffusiveFlux = claw.newDiffusiveFlux(cacheMap);
// compute the gradient of solution if necessary
fullMatrix<double> *gradientSolutionQP= 0;
if (gradSolutionQPe.somethingDependOnMe()) {
gradientSolutionQP = new fullMatrix<double> (group.getNbIntegrationPoints(), group.getNbElements() * nbFields * 3);
//group.getGradientOfSolution().mult ( group.getCollocationMatrix() , *gradientSolutionQP);
}
fullMatrix<double> fuvwe;
fullMatrix<double> source;
fullMatrix<double> dPsiDx,dofs;
// ----- 2.3 --- iterate on elements
for (int iElement=0 ; iElement<group.getNbElements() ;++iElement) {
// ----- 2.3.1 --- build a small object that contains elementary solution, jacobians, gmsh element
solutionQPe.setAsProxy(solutionQP, iElement*nbFields, nbFields );
if (gradSolutionQPe.somethingDependOnMe())
gradSolutionQPe.setAsProxy(*gradientSolutionQP, 3*iElement*nbFields,3*nbFields );
if(gradientSolutionQPe.somethingDependOnMe()){
dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
dofs.setAsProxy(solution, nbFields*iElement, nbFields);
dPsiDx.mult(dofs, gradientSolutionQPe.set());
}
cacheElement.set(group.getElement(iElement));
if(convectiveFlux || diffusiveFlux) {
// ----- 2.3.3 --- compute fluxes in UVW coordinates
......@@ -84,13 +83,14 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
for (int k=0;k<nbFields;k++){
if(convectiveFlux)
fuvwe(iPt,k) += ( (*convectiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
if(diffusiveFlux)
if(diffusiveFlux){
fuvwe(iPt,k) += ( (*diffusiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
}
}
}
}
}
}
if (sourceTerm){
source.setAsProxy(Source, iElement*claw.nbFields(),claw.nbFields());
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
......@@ -100,8 +100,6 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
}
}
if(gradientSolutionQP)
delete gradientSolutionQP;
}
// ----- 3 ---- do the redistribution at nodes using as many BLAS3 operations as there are local coordinates
if(convectiveFlux || diffusiveFlux){
......@@ -158,8 +156,8 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
dataCacheDouble &uvwRight = cacheMapRight.provideData("UVW");
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
dataCacheDouble &solutionQPRight = cacheMapRight.provideData("Solution");
dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("GradientSolution");
dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("GradientSolution");
dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient");
dataCacheDouble &gradientSolutionRight = cacheMapRight.provideData("SolutionGradient");
//resize gradientSolutionLeft and Right
gradientSolutionLeft.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
gradientSolutionRight.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
......@@ -167,9 +165,14 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
// A2 ) ask the law to provide the fluxes (possibly NULL)
dataCacheDouble *riemannSolver = claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
dataCacheDouble *maximumDiffusivityLeft = claw.newMaximumDiffusivity(cacheMapLeft);
dataCacheDouble *maximumDiffusivityRight = claw.newMaximumDiffusivity(cacheMapRight);
dataCacheDouble *diffusiveFluxRight = claw.newDiffusiveFlux(cacheMapRight);
fullMatrix<double> dPsiLeftDx,dPsiRightDx,dofsLeft,dofsRight,normalFluxQP;
int p = group.getGroupLeft().getOrder();
int dim = group.getGroupLeft().getElement(0)->getDim();
for (int iFace=0 ; iFace < nbFaces ; ++iFace) {
// B1 ) adjust the proxies for this element
// NB : the B1 section will almost completely disapear
......@@ -183,7 +186,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
normals.setAsProxy(group.getNormals(), iFace*group.getNbIntegrationPoints(),group.getNbIntegrationPoints());
// proxies needed to compute the gradient of the solution and the IP term
dPsiLeftDx.setAsProxy(DPsiLeftDx,iFace*nbNodesLeft,nbNodesLeft);
dPsiRightDx.setAsProxy(DPsiLeftDx,iFace*nbNodesRight,nbNodesRight);
dPsiRightDx.setAsProxy(DPsiRightDx,iFace*nbNodesRight,nbNodesRight);
dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
dofsRight.setAsProxy(solutionRight, nbFields*group.getElementRightId(iFace), nbFields);
uvwLeft.setAsProxy(group.getIntegrationOnElementLeft(iFace));
......@@ -198,27 +201,44 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
}
// B3 ) compute fluxes
//JF : here you can test (diffusiveFluxLeft!=NULL) to know if the law is diffusive
//you can access the values by (*diffusiveFluxLeft)(iQP*3+iXYZ, iField)
//use gradientSolutionLef(IQP*3+IXYZ, iField) and dPsiLeftDx(iQP*3+iXYZ, iPsi)
if (diffusiveFluxLeft) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
//just for the lisibility :
const fullMatrix<double> &dfl = (*diffusiveFluxLeft)();
const fullMatrix<double> &dfr = (*diffusiveFluxRight)();
for (int k=0;k<nbFields;k++) {
double meanNormalFlux =
((dfl(iPt,k*3 )+dfr(iPt,k*3 )) * normals(0,iPt)
+(dfl(iPt,k*3+1)+dfr(iPt,k*3+1)) * normals(1,iPt)
+(dfl(iPt,k*3+2)+dfr(iPt,k*3+2)) * normals(2,iPt))/2;
double minl = std::min(group.getElementVolumeLeft(iFace),
group.getElementVolumeRight(iFace)
)/group.getInterfaceSurface(iFace);
double nu = std::max((*maximumDiffusivityRight)(iPt,0),(*maximumDiffusivityLeft)(iPt,0));
double mu = (p+1)*(p+dim)/dim*nu/minl;
double solutionJumpPenalty = (solutionQPLeft(iPt,k)-solutionQPRight(iPt,k))/2*mu;
normalFluxQP(iPt,k) -= (meanNormalFlux+solutionJumpPenalty)*detJ;
normalFluxQP(iPt,k+nbFields) += (meanNormalFlux+solutionJumpPenalty)*detJ;
}
}
}
if (riemannSolver) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<nbFields*2;k++)
normalFluxQP(iPt,k) = (*riemannSolver)(iPt,k)*detJ;
normalFluxQP(iPt,k) += (*riemannSolver)(iPt,k)*detJ;
}
}
}
// C ) redistribute the flux to the residual (at Faces nodes)
if(riemannSolver || diffusiveFluxLeft)
residual.gemm(group.getRedistributionMatrix(),NormalFluxQP);
// D ) delete the dataCacheDouble provided by the law
if(riemannSolver)
delete riemannSolver;
if(diffusiveFluxLeft) {
delete diffusiveFluxLeft;
delete diffusiveFluxRight;
}
}
void dgAlgorithm::multAddInverseMassMatrix ( /*dofManager &dof,*/
......
......@@ -15,7 +15,7 @@ class dgAlgorithm {
static void residualVolume ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
const dgGroupOfElements & group,
const fullMatrix<double> &solution,
fullMatrix<double> &solution,
fullMatrix<double> &residual);
static void residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
const dgConservationLaw &claw, // the conservation law
......
......@@ -32,6 +32,8 @@ class dgConservationLaw {
virtual dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newDiffusiveFlux (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newDiffusiveInterfaceTerm (dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {return NULL;}
virtual dataCacheDouble *newMaximumDiffusivity (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newConvectiveFlux (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const {return NULL;}
virtual dataCacheDouble *newRiemannSolver (dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {return NULL;}
......@@ -51,7 +53,7 @@ class dgConservationLaw {
};
dgConservationLaw *dgNewConservationLawAdvection(const std::string vname);
dgConservationLaw *dgNewConservationLawAdvection(const std::string vname,const std::string nuname);
dgConservationLaw *dgNewConservationLawShallowWater2d();
dgConservationLaw *dgNewConservationLawWaveEquation();
dgConservationLaw *dgNewPerfectGasLaw2d();
......
......@@ -6,7 +6,7 @@
#include "function.h"
class dgConservationLawAdvection : public dgConservationLaw {
std::string _vFunctionName;
std::string _vFunctionName,_nuFunctionName;
class advection : public dataCacheDouble {
dataCacheDouble &sol, &v;
public:
......@@ -46,23 +46,56 @@ class dgConservationLawAdvection : public dgConservationLaw {
}
}
};
class diffusion : public dataCacheDouble {
dataCacheDouble &solgrad, &nu;
public:
diffusion(std::string nuFunctionName, dataCacheMap &cacheMap):
solgrad(cacheMap.get("SolutionGradient",this)),
nu(cacheMap.get(nuFunctionName,this))
{};
void _eval () {
if(_value.size1() != solgrad().size1()/3)
_value=fullMatrix<double>(solgrad().size1()/3,3);
for(int i=0; i< solgrad().size1()/3; i++) {
_value(i,0) = -solgrad(i*3,0)*nu(i,0);
_value(i,1) = -solgrad(i*3+1,0)*nu(i,0);
_value(i,2) = -solgrad(i*3+1,0)*nu(i,0);
}
}
};
public:
dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const {
if( !_vFunctionName.empty())
return new advection(_vFunctionName,cacheMap);
else
return NULL;
}
dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const {
if( !_nuFunctionName.empty())
return &cacheMap.get(_nuFunctionName);
else
return NULL;
}
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const {
if( !_vFunctionName.empty())
return new riemann(_vFunctionName,cacheMapLeft, cacheMapRight);
else
return NULL;
}
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const {
return 0;
if( !_nuFunctionName.empty())
return new diffusion(_nuFunctionName,cacheMap);
else
return NULL;
}
dgConservationLawAdvection(std::string vFunctionName)
dgConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName)
{
_vFunctionName = vFunctionName;
_nuFunctionName = nuFunctionName;
_nbf = 1;
}
};
dgConservationLaw *dgNewConservationLawAdvection(std::string vFunctionName) {
return new dgConservationLawAdvection(vFunctionName);
dgConservationLaw *dgNewConservationLawAdvection(std::string vFunctionName, std::string nuFunctionName) {
return new dgConservationLawAdvection(vFunctionName,nuFunctionName);
}
......@@ -47,7 +47,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
_fs(*_elements[0]->getFunctionSpace(polyOrder)),
_integration(dgGetIntegrationRule (_elements[0], polyOrder))
{
_order=polyOrder;
_dimUVW = _dimXYZ = e[0]->getDim();
// this is the biggest piece of data ... the mappings
int nbPsi = _fs.coefficients.size1();
......@@ -59,6 +59,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
_mapping = new fullMatrix<double> (e.size(), 10 * _integration->size1());
_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);
double g[256][3],f[256];
for (int i=0;i<_elements.size();i++){
......@@ -70,6 +71,7 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
double jac[3][3],ijac[3][3],detjac;
(*_mapping)(i,10*j + 9) =
e->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
(*_elementVolume)(i,0) += (*_mapping)(i,10*j+9)*(*_integration)(j,3);
const double weight = (*_integration)(j,3);
detjac=inv3x3(jac,ijac);
(*_mapping)(i,10*j + 0) = ijac[0][0];
......@@ -123,6 +125,7 @@ dgGroupOfElements::~dgGroupOfElements(){
delete _mapping;
delete _collocation;
delete _imass;
delete _elementVolume;
}
......@@ -204,6 +207,7 @@ void dgGroupOfFaces::init(int pOrder) {
_redistribution = new fullMatrix<double> (_fsFace->coefficients.size1(),_integration->size1());
_collocation = new fullMatrix<double> (_integration->size1(), _fsFace->coefficients.size1());
_detJac = new fullMatrix<double> (_integration->size1(), _faces.size());
_interfaceSurface = new fullMatrix<double>(_faces.size(),1);
for (size_t i=0; i<_closuresLeft.size(); i++)
_integrationPointsLeft.push_back(dgGetFaceIntegrationRuleOnElement(_fsFace,*_integration,_fsLeft,_closuresLeft[i]));
for (size_t i=0; i<_closuresRight.size(); i++)
......@@ -222,6 +226,7 @@ void dgGroupOfFaces::init(int pOrder) {
for (int j=0;j< _integration->size1() ; j++ ){
double jac[3][3];
(*_detJac)(j,i) = f->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
(*_interfaceSurface)(i,0) += (*_integration)(j,3)*(*_detJac)(j,i);
}
}
......@@ -296,6 +301,7 @@ dgGroupOfFaces::~dgGroupOfFaces()
delete _normals;
delete _dPsiLeftDxOnQP;
delete _dPsiRightDxOnQP;
delete _interfaceSurface;
}
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices):
......
......@@ -63,6 +63,11 @@ class dgGroupOfElements {
// dimension of the parametric space and of the real space
// may be different if the domain is a surface in 3D (manifold)
int _dimUVW, _dimXYZ;
// polynomial order of the interpolation
int _order;
//
// volume/surface/length of the element (sum_qp w_i detJ_i)
fullMatrix<double> *_elementVolume;
// forbid the copy
// dgGroupOfElements (const dgGroupOfElements &e, int order) {}
// dgGroupOfElements & operator = (const dgGroupOfElements &e) {}
......@@ -79,11 +84,13 @@ public:
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
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 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;}
inline fullMatrix<double> &getInverseMassMatrix () const {return *_imass;}
inline const fullMatrix<double> getMapping (int iElement) const {return fullMatrix<double>(*_mapping, iElement, 1);}
inline int getOrder() const {return _order;}
};
class dgGroupOfFaces {
......@@ -127,6 +134,8 @@ class dgGroupOfFaces {
//fullMatrix<double> *_collocationLeft, *_collocationRight;
// redistribution matrices \psi_i (GP_j) * weight_j
fullMatrix<double> *_redistribution;
// surface/length/1 of the interface element (sum_qp w_i detJ_i)
fullMatrix<double> *_interfaceSurface;
//common part of the 3 constructors
void init(int pOrder);
public:
......@@ -136,6 +145,8 @@ public:
inline int getElementRightId (int i) const {return _right[i];};
inline MElement* getElementLeft (int i) const {return _groupLeft.getElement(_left[i]);}
inline MElement* getElementRight (int i) const {return _groupRight.getElement(_right[i]);}
inline double getElementVolumeLeft(int iFace) const {return _groupLeft.getElementVolume(_left[iFace]);}
inline double getElementVolumeRight(int iFace) const {return _groupRight.getElementVolume(_right[iFace]);}
inline MElement* getFace (int iElement) const {return _faces[iElement];}
inline const std::vector<int> &getClosureLeft(int iFace) const{ return _closuresLeft[_closuresIdLeft[iFace]];}
inline const std::vector<int> &getClosureRight(int iFace) const{ return _closuresRight[_closuresIdRight[iFace]];}
......@@ -160,6 +171,7 @@ public:
inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getRedistributionMatrix () const {return *_redistribution;}
inline double getDetJ (int iElement, int iGaussPoint) const {return (*_detJac)(iGaussPoint,iElement);}
inline double getInterfaceSurface (int iFace)const {return (*_interfaceSurface)(iFace,0);}
//keep this outside the Algorithm because this is the only place where data overlap
void mapToInterface(int nFields, const fullMatrix<double> &vLeft, const fullMatrix<double> &vRight, fullMatrix<double> &v);
void mapFromInterface(int nFields, const fullMatrix<double> &v, fullMatrix<double> &vLeft, fullMatrix<double> &vRight);
......
......@@ -69,8 +69,7 @@ int dgSystemOfEquations::setConservationLaw(lua_State *L){
else if (_cLawName == "PerfectGas2d")
_claw = dgNewPerfectGasLaw2d();
else if (_cLawName == "AdvectionDiffusion"){
std::string advFunction = luaL_checkstring(L,2);
_claw = dgNewConservationLawAdvection(advFunction);
_claw = dgNewConservationLawAdvection(luaL_checkstring(L,2),luaL_checkstring(L,3));
}
if (!_claw)throw;
return 0;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment