diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 5268a3148653124b18f7145947c62dd075efe57b..60eb285d8e0e556db72ddb30381bc8edf1ab75d4 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -59,6 +59,7 @@ void PrintUsage(const char *name)
Msg::Direct(" -1, -2, -3 Perform 1D, 2D or 3D mesh generation, then exit");
Msg::Direct(" -refine Perform uniform mesh refinement, then exit");
Msg::Direct(" -part int Partition after batch mesh generation");
+ Msg::Direct(" -renumber Renumber the mesh elements after batch mesh generation");
Msg::Direct(" -saveall Save all elements (discard physical group definitions)");
Msg::Direct(" -o file Specify mesh output file name");
Msg::Direct(" -format string Set output mesh format (msh, msh1, msh2, unv, vrml, stl, mesh,");
@@ -170,10 +171,15 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->batch = 5;
i++;
}
+ else if(!strcmp(argv[i] + 1, "renumber")) {
+ CTX::instance()->batchAfterMesh = 1;
+ CTX::instance()->partitionOptions.renumber = 1;
+ i++;
+ }
else if(!strcmp(argv[i] + 1, "part")) {
i++;
if(argv[i]){
- CTX::instance()->batchAfterMesh = 1;
+ CTX::instance()->batchAfterMesh =1 ;
opt_mesh_partition_num(0, GMSH_SET, atoi(argv[i++]));
}
else
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index ae9e4b86c40ed71cbdae59eb6978ee2b2b1ac84b..fca933aaf72b6ab082fc8680f25fdf9860688d40 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1109,7 +1109,7 @@ StringXNumber MeshOptions_Number[] = {
"Number of hexahedra in the current mesh (read-only)" },
{ F, "NbNodes" , opt_mesh_nb_nodes , 0. ,
"Number of nodes in the current mesh (read-only)" },
- { F|O, "NbPartitions" , opt_mesh_partition_num, 4. ,
+ { F|O, "NbPartitions" , opt_mesh_partition_num, 1. ,
"Number of partitions" },
{ F, "NbPrisms" , opt_mesh_nb_prisms , 0. ,
"Number of prisms in the current mesh (read-only)" },
diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index cbcb48aa248e40e56061e086b7361d30c50a0c82..25d5fe73b46005fc0ef35ba886e66a4a987abd1e 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -169,8 +169,10 @@ int GmshBatch()
else if(CTX::instance()->batch == 5)
RefineMesh(GModel::current(), CTX::instance()->mesh.secondOrderLinear);
#if defined(HAVE_CHACO) || defined(HAVE_METIS)
- if(CTX::instance()->batchAfterMesh == 1)
- PartitionMesh(GModel::current(), CTX::instance()->partitionOptions);
+ if(CTX::instance()->batchAfterMesh == 1){
+ if (CTX::instance()->partitionOptions.num_partitions > 1)PartitionMesh(GModel::current(), CTX::instance()->partitionOptions);
+ if (CTX::instance()->partitionOptions.renumber)RenumberMesh(GModel::current(), CTX::instance()->partitionOptions);
+ }
#endif
#endif
CreateOutputFile(CTX::instance()->outputFileName, CTX::instance()->mesh.format);
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 4d2f3fde2ee2040ee621a9939f4c8dca633d2723..25f0f3cc864a4ef75b6ffd56d80b824e6cbbd731 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1234,7 +1234,7 @@ void deMeshGFace::operator() (GFace *gf)
gf->meshStatistics.nbTriangle = gf->meshStatistics.nbEdge = 0;
}
-const int debugSurface = -1;
+int debugSurface = -1;
void meshGFace::operator() (GFace *gf)
{
diff --git a/Mesh/meshPartition.cpp b/Mesh/meshPartition.cpp
index 3b9cf5e2385764cf721d3ef234d1baa709a22da9..246839dbe3b98760407b0fcb9d14488795880ab0 100644
--- a/Mesh/meshPartition.cpp
+++ b/Mesh/meshPartition.cpp
@@ -52,6 +52,7 @@ extern "C" void METIS_PartGraphRecursive
extern "C" void METIS_PartGraphKway
(int *, idxtype *, idxtype *, idxtype *, idxtype *, int *, int *, int *, int *,
int *, idxtype *);
+extern "C" void METIS_NodeND(int *n,int *xadj,int *adjncy,int *numflag,int *options,int *perm,int *iperm);
/*==============================================================================
@@ -112,6 +113,34 @@ void MakeGraphDIM(const EntIter begin, const EntIter end,
+int RenumberMesh(GModel *const model, meshPartitionOptions &options, std::vector<MElement*> &numbered)
+{
+
+ Graph graph;
+ BoElemGrVec boElemGrVec;
+ int ier;
+ Msg::StatusBar(1, true, "Building graph...");
+ ier = MakeGraph(model, graph, &boElemGrVec);
+ Msg::StatusBar(1, true, "Renumbering graph...");
+ if(!ier) ier = RenumberGraph(graph, options);
+ if(ier) {
+ Msg::StatusBar(1, false, "Mesh");
+ return 1;
+ }
+
+ // create the numbering
+ numbered.clear();
+ const int n = graph.getNumVertex();
+ numbered.resize(n);
+ for(int i = 0; i != n; ++i) {
+ numbered[graph.partition[i]-1] = graph.element[i];
+ }
+
+ Msg::Info("Renumbering complete");
+ Msg::StatusBar(1, false, "Mesh");
+ return 0;
+}
+
int PartitionMeshElements(std::vector<MElement*> &elements, meshPartitionOptions &options)
{
GModel *tmp_model = new GModel();
@@ -153,6 +182,80 @@ int PartitionMeshFace(std::list<GFace*> &cFaces, meshPartitionOptions &options)
delete tmp_model;
return 1;
}
+int RenumberMeshElements( std::vector<MElement*> &elements, meshPartitionOptions &options ){
+ if (elements.size() < 3)return 1;
+ GModel *tmp_model = new GModel();
+ std::set<MVertex *> setv;
+ for (unsigned i=0;i<elements.size();++i)
+ for (int j=0;j<elements[i]->getNumVertices();j++)
+ setv.insert(elements[i]->getVertex(j));
+ if (elements[0]->getDim() == 2){
+ GFace *gf = new discreteFace(tmp_model, 1);
+ for (std::set<MVertex* >::iterator it = setv.begin(); it != setv.end(); it++)
+ gf->mesh_vertices.push_back(*it);
+ for (std::vector<MElement* >::iterator it = elements.begin(); it != elements.end(); it++){
+ if ((*it)->getType() == TYPE_TRI)
+ gf->triangles.push_back((MTriangle*)(*it));
+ else if ((*it)->getType() == TYPE_QUA)
+ gf->quadrangles.push_back((MQuadrangle*)(*it));
+ }
+ tmp_model->add(gf);
+ RenumberMesh(tmp_model,options,elements);
+ tmp_model->remove(gf);
+ }
+ else if (elements[0]->getDim() == 3){
+ GFace *gf = new discreteFace(tmp_model, 1);
+ for (std::set<MVertex* >::iterator it = setv.begin(); it != setv.end(); it++)
+ gf->mesh_vertices.push_back(*it);
+ for (std::vector<MElement* >::iterator it = elements.begin(); it != elements.end(); it++){
+ if ((*it)->getType() == TYPE_TRI)
+ gf->triangles.push_back((MTriangle*)(*it));
+ else if ((*it)->getType() == TYPE_QUA)
+ gf->quadrangles.push_back((MQuadrangle*)(*it));
+ }
+ tmp_model->add(gf);
+ }
+ delete tmp_model;
+ return 1;
+}
+
+int RenumberMesh(GModel *const model, meshPartitionOptions &options)
+{
+
+ for (GModel::fiter it = model->firstFace() ; it != model->lastFace() ; ++it){
+ std::vector<MElement *> temp;
+
+ temp.insert(temp.begin(),(*it)->triangles.begin(),(*it)->triangles.end());
+ RenumberMeshElements (temp,options);
+ (*it)->triangles.clear();
+ for (int i=0;i<temp.size();i++) (*it)->triangles.push_back((MTriangle*)temp[i]);
+
+ temp.clear();
+
+ temp.insert(temp.begin(),(*it)->quadrangles.begin(),(*it)->quadrangles.end());
+ RenumberMeshElements (temp,options);
+ (*it)->quadrangles.clear();
+ for (int i=0;i<temp.size();i++) (*it)->quadrangles.push_back((MQuadrangle*)temp[i]);
+ }
+ for (GModel::riter it = model->firstRegion() ; it != model->lastRegion() ; ++it){
+ std::vector<MElement *> temp;
+
+ temp.insert(temp.begin(),(*it)->tetrahedra.begin(),(*it)->tetrahedra.end());
+ RenumberMeshElements (temp,options);
+ (*it)->tetrahedra.clear();
+ for (int i=0;i<temp.size();i++) (*it)->tetrahedra.push_back((MTetrahedron*)temp[i]);
+
+ temp.clear();
+
+ temp.insert(temp.begin(),(*it)->hexahedra.begin(),(*it)->hexahedra.end());
+ RenumberMeshElements (temp,options);
+ (*it)->hexahedra.clear();
+ for (int i=0;i<temp.size();i++) (*it)->hexahedra.push_back((MHexahedron*)temp[i]);
+ }
+ return 1;
+}
+
+
int PartitionMesh(GModel *const model, meshPartitionOptions &options)
{
@@ -325,7 +428,7 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
&options.num_partitions, metisOptions, &edgeCut,
&graph.partition[graph.section[iSec]]);
break;
- }
+ }
}
catch(...) {
Msg::Error("METIS threw an exception");
@@ -343,6 +446,41 @@ int PartitionGraph(Graph &graph, meshPartitionOptions &options)
return ier;
}
+int RenumberGraph(Graph &graph, meshPartitionOptions &options)
+{
+ int ier = 0;
+#ifdef HAVE_METIS
+ {
+ Msg::Info("Launching METIS graph renumberer");
+ // "C" numbering for Metis
+ {
+ int *p = &graph.adjncy[0]; //**Sections
+ for(int n = graph.adjncy.size(); n--;) --(*p++);
+ }
+ try {
+ int n = graph.getNumVertex();
+ int numflag = 0;
+ const int iSec = 0;
+ int options = 0;
+ int *perm = new int[n];
+ METIS_NodeND(&n,
+ &graph.xadj[graph.section[iSec]],
+ &graph.adjncy[graph.section[iSec]], &numflag,&options,perm,&graph.partition[graph.section[iSec]]);
+ delete [] perm;
+ }
+ catch(...) {
+ Msg::Error("METIS threw an exception");
+ ier = 2;
+ }
+ // Number elements from 1
+ if(!ier) {
+ int *p = &graph.partition[0]; //**Sections
+ for(int n = graph.getNumVertex(); n--;) ++(*p++);
+ }
+ }
+#endif
+ return ier;
+}
/*******************************************************************************
*
diff --git a/Mesh/meshPartition.h b/Mesh/meshPartition.h
index ac97f350f873233139c7702a1727b302a32ce2ad..240e9e4b50330e8b7523f97b97922875cf86b5b6 100644
--- a/Mesh/meshPartition.h
+++ b/Mesh/meshPartition.h
@@ -23,7 +23,9 @@ struct meshPartitionOptions;
******************************************************************************/
int PartitionGraph(Graph &graph, meshPartitionOptions &options);
+int RenumberGraph(Graph &graph, meshPartitionOptions &options);
int PartitionMesh(GModel *const model, meshPartitionOptions &options);
+int RenumberMesh(GModel *const model, meshPartitionOptions &options);
int PartitionMeshFace(std::list<GFace*> &cFaces, meshPartitionOptions &options);
int PartitionMeshElements(std::vector<MElement*> &elements, meshPartitionOptions &options);
bool PartitionZeroGenus(std::list<GFace*> &cFaces, int &nbParts);
diff --git a/Mesh/meshPartitionOptions.h b/Mesh/meshPartitionOptions.h
index 4488a8baa5387c4bd8525a8715bc11e102987e1b..05406a7715b44f5f4d547cdd2a0e7558da7c29de 100644
--- a/Mesh/meshPartitionOptions.h
+++ b/Mesh/meshPartitionOptions.h
@@ -14,6 +14,8 @@ struct meshPartitionOptions
// 2 - METIS
int num_partitions;
+ int renumber;
+
bool createPartitionBoundaries;
//--Chaco
@@ -82,7 +84,8 @@ struct meshPartitionOptions
void setDefaults()
{
partitioner = 2;
- num_partitions = 4;
+ num_partitions = 1;
+ renumber = 0;
global_method = 1;
architecture = 1;
ndims_tot = 2;
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index c17881309e0edc37854c0eff365569d9ccd1aaca..dcbb2416a9aa5230eed6b3b68813879020b77d68 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -88,6 +88,12 @@ class fullVector
}
printf("\n");
}
+ void binarySave (FILE *f) const{
+ fwrite (_data,sizeof(scalar),_r,f);
+ }
+ void binaryLoad (FILE *f){
+ fread (_data,sizeof(scalar),_r,f);
+ }
};
template <class scalar>
diff --git a/Solver/TESTCASES/CylinderEddies.lua b/Solver/TESTCASES/CylinderEddies.lua
new file mode 100644
index 0000000000000000000000000000000000000000..cd73f5cf70b3b8a22126b1ed28576d79de2aa293
--- /dev/null
+++ b/Solver/TESTCASES/CylinderEddies.lua
@@ -0,0 +1,79 @@
+MACH = .05;
+RHO = 1.0;
+PRES = 1./(MACH*RHO*RHO*1.4*1.4)
+V = 1.0
+SOUND = V/MACH
+
+--[[
+ Function for initial conditions
+--]]
+function free_stream( XYZ, FCT )
+ for i=0,XYZ:size1()-1 do
+ FCT:set(i,0,RHO)
+ FCT:set(i,1,RHO*V)
+ FCT:set(i,2,0.0)
+ FCT:set(i,3, 0.5*RHO*V*V+PRES/0.4)
+ end
+end
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+order = 3
+
+FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
+
+-- diffusivity
+mu=fullMatrix(1,1);
+mu:set(0,0,0.03)
+kappa=fullMatrix(1,1);
+kappa:set(0,0,0.03)
+
+print'*** Loading the mesh and the model ***'
+myModel = GModel ()
+myModel:load ('cyl.geo')
+myModel:load ('cyl.msh')
+
+print'*** Create a dg solver ***'
+DG = dgSystemOfEquations (myModel)
+DG:setOrder(order)
+law=dgPerfectGasLaw2d()
+
+law:setViscosityAndThermalConductivity(functionConstant(mu):getName(),functionConstant(kappa):getName());
+
+DG:setConservationLaw(law)
+law:addBoundaryCondition('Cylinder',law:newNonSlipWallBoundary())
+law:addBoundaryCondition('Box' ,law:newOutsideValueBoundary(FS))
+
+DG:setup()
+
+print'*** setting the initial solution ***'
+
+DG:L2Projection(FS)
+
+--print'*** export ***'
+
+DG:exportSolution('output/cyl_0')
+
+print'*** solve ***'
+
+CFL = 2.1;
+dt = CFL * DG:computeInvSpectralRadius();
+print('DT = ',dt)
+T = 0;
+for i=1,10000 do
+ dt = CFL * DG:computeInvSpectralRadius();
+ norm = DG:RK44(dt)
+ T = T + dt
+ if (i % 1 == 0) then
+ print('*** ITER ***',i,norm,dt,T)
+ end
+ if (i % 100 == 0) then
+ DG:exportSolution(string.format("output/cyl-%06d", i))
+ end
+end
+
+print'*** done ***'
+
+
diff --git a/Solver/TESTCASES/RayleighTaylor.lua b/Solver/TESTCASES/RayleighTaylor.lua
new file mode 100644
index 0000000000000000000000000000000000000000..f58f88d7762b500864b89a6165e864a1292a835b
--- /dev/null
+++ b/Solver/TESTCASES/RayleighTaylor.lua
@@ -0,0 +1,103 @@
+ATWOOD = .5
+RHO1 = 1.0
+RHO2 = RHO1 * (1.+ATWOOD)/(1.-ATWOOD)
+taper = 6.0
+--Mach number of perturbation
+machnum =0.1
+GAMMA = 1.4
+-- Min speed of sound is at top where p=1, rho=2
+epsilon_z= machnum * math.sqrt(GAMMA*1.0/2.0)
+epsilon_xy= -epsilon_z*taper/16.0
+BANDWIDTH = 0.005;
+PI = 3.14159
+
+--[[
+ Function for initial conditions
+--]]
+function initial_condition( XYZ, FCT )
+ for i=0,XYZ:size1()-1 do
+ X = XYZ:get(i,0)-.125
+ Y = XYZ:get(i,1)
+ RHO = 1;
+ if Y > 0.0 then
+ RHO = 2
+ end
+-- smooth version
+ RHO = 1.0 + (math.atan (Y/BANDWIDTH) / PI + 0.5)
+
+ PRES = 2 - RHO*Y;
+ vx = epsilon_xy * math.sin(X/0.25 * 2*PI) * math.cos(Y*PI) * math.pow(math.cos(Y*PI),taper-1.0);
+ vy = epsilon_z * math.cos(X/0.25 * 2*PI) * math.pow(math.cos(Y*PI),taper);
+ FCT:set(i,0,RHO)
+ FCT:set(i,1,RHO*vx)
+ FCT:set(i,2,RHO*vy)
+ FCT:set(i,3, 0.5*RHO*(vx*vx+vy*vy)+PRES/0.4)
+ end
+end
+
+--[[
+ Example of a lua program driving the DG code
+--]]
+
+order = 4
+print'*** Loading the mesh and the model ***'
+myModel = GModel ()
+ myModel:load ('rect.geo')
+if (order == 1) then
+ myModel:load ('rect.msh')
+elseif (order == 2) then
+ myModel:load ('rect2.msh')
+elseif (order == 3) then
+ myModel:load ('rect2.msh')
+elseif (order == 4) then
+ myModel:load ('rect4.msh')
+end
+
+print'*** Create a dg solver ***'
+DG = dgSystemOfEquations (myModel)
+DG:setOrder(order)
+law=dgPerfectGasLaw2d()
+
+g=fullMatrix(4,1);
+g:set(0,0,0)
+g:set(1,0,0)
+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)
+DG:limitSolution()
+
+--print'*** export ***'
+
+DG:exportSolution('output/rt_0')
+
+print'*** solve ***'
+
+LC = 0.1*.1
+dt = .0003;
+print('DT=',dt)
+
+for i=1,10000 do
+-- norm = DG:RK44_limiter(dt)
+ norm = DG:RK44(dt)
+ if (i % 10 == 0) then
+ print('*** ITER ***',i,norm)
+ end
+ if (i % 100 == 0) then
+ DG:exportSolution(string.format("output/rt-%06d", i))
+ end
+end
+
+print'*** done ***'
+
+
diff --git a/Solver/TESTCASES/cyl.geo b/Solver/TESTCASES/cyl.geo
new file mode 100644
index 0000000000000000000000000000000000000000..efd6f094efcc21fcfb8cdd95d2e2d96b2bcd9706
--- /dev/null
+++ b/Solver/TESTCASES/cyl.geo
@@ -0,0 +1,25 @@
+a = .1;
+b = .5;
+Point(1) = {0, 0, 0, a};
+Point(2) = {1, 0, 0, a};
+Point(3) = {-1, 0, 0, a};
+Point(4) = {0, 1, 0, a};
+Point(5) = {0, -1, 0, a};
+Point(6) = {-5, -5, 0, b};
+Point(7) = {-5, 5, 0, b};
+Point(8) = {15, 5, 0, b};
+Point(9) = {15, -5, 0, b};
+Circle(1) = {4, 1, 2};
+Circle(2) = {2, 1, 5};
+Circle(3) = {5, 1, 3};
+Circle(4) = {3, 1, 4};
+Line(5) = {7, 8};
+Line(6) = {8, 9};
+Line(7) = {9, 6};
+Line(8) = {6, 7};
+Line Loop(9) = {6, 7, 8, 5};
+Line Loop(10) = {1, 2, 3, 4};
+Plane Surface(11) = {9, 10};
+Physical Surface("Air") = {11};
+Physical Line("Cylinder") = {4, 1, 2, 3};
+Physical Line("Box") = {5, 6, 7, 8};
diff --git a/Solver/TESTCASES/rect.geo b/Solver/TESTCASES/rect.geo
new file mode 100644
index 0000000000000000000000000000000000000000..7deb863a25bf03481ba53710aa435a9b9be0b603
--- /dev/null
+++ b/Solver/TESTCASES/rect.geo
@@ -0,0 +1,25 @@
+Point(1) = {-.125, 0, 0, .05};
+Point(2) = {.125, 0, 0, .05};
+Point(3) = {.125, .5, 0, .05};
+Point(4) = {-.125, .5, 0, .05};
+Point(5) = {-.125, -.5, 0, .05};
+Point(6) = {.125, -.5, 0, .05};
+Line(1) = {5, 1};
+Line(2) = {1, 2};
+Line(3) = {2, 6};
+Line(4) = {6, 5};
+Line(5) = {1, 4};
+Line(6) = {4, 3};
+Line(7) = {3, 2};
+Line Loop(8) = {7, -2, 5, 6};
+Plane Surface(9) = {8};
+Line Loop(10) = {3, 4, 1, 2};
+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 Surface {9} Alternated;
+Transfinite Surface {11} Alternated;
+Recombine Surface {9, 11};
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index a86a0d4305f7da8090feaf1670215968bcee8253..24a9de318ef0394b1b952c03ca66c3fd8403f57d 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -9,6 +9,7 @@
#include "MEdge.h"
#include "function.h"
#include "dgLimiter.h"
+#include "meshPartition.h"
/*
compute
@@ -62,6 +63,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
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(gradientSolutionQPe.somethingDependOnMe()){
dPsiDx.setAsProxy(group.getDPsiDx(),iElement*group.getNbNodes(),group.getNbNodes());
dofs.setAsProxy(solution, nbFields*iElement, nbFields);
@@ -86,7 +88,7 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
if(convectiveFlux)
fuvwe(iPt,k) += ( (*convectiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
if(diffusiveFlux){
- fuvwe(iPt,k) += ( (*diffusiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
+ fuvwe(iPt,k) += ( (*diffusiveFlux)(iPt,iXYZ*nbFields+k)) * factor;
}
}
}
@@ -202,7 +204,7 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
dPsiRightDx.mult(dofsRight, gradientSolutionRight.set());
}
- // B3 ) compute fluxes
+ // B3 ) compute fluxes
if (diffusiveFluxLeft) {
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
@@ -211,9 +213,9 @@ void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (may
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;
+ ((dfl(iPt,k+nbFields*0 )+dfr(iPt,k+nbFields*0)) * normals(0,iPt)
+ +(dfl(iPt,k+nbFields*1 )+dfr(iPt,k+nbFields*1)) * normals(1,iPt)
+ +(dfl(iPt,k+nbFields*2 )+dfr(iPt,k+nbFields*2)) * normals(2,iPt))/2;
double minl = std::min(group.getElementVolumeLeft(iFace),
group.getElementVolumeRight(iFace)
)/group.getInterfaceSurface(iFace);
@@ -467,6 +469,7 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
)
{
int nbFields = claw.nbFields();
+ int nbNodesLeft = group.getGroupLeft().getNbNodes();
const dgBoundaryCondition *boundaryCondition = claw.getBoundaryCondition(group.getBoundaryTag());
// ----- 1 ---- get the solution at quadrature points
fullMatrix<double> solutionQP (group.getNbIntegrationPoints(),group.getNbElements() * nbFields);
@@ -475,33 +478,73 @@ void dgAlgorithm::residualBoundary ( //dofManager &dof, // the DOF manager (mayb
fullMatrix<double> NormalFluxQP ( group.getNbIntegrationPoints(), group.getNbElements()*nbFields);
dataCacheMap cacheMapLeft(group.getNbIntegrationPoints());
+ fullMatrix<double> &DPsiLeftDx = group.getDPsiLeftDxMatrix();
// provided dataCache
dataCacheDouble &uvw=cacheMapLeft.provideData("UVW");
dataCacheDouble &solutionQPLeft = cacheMapLeft.provideData("Solution");
+ dataCacheDouble &gradientSolutionLeft = cacheMapLeft.provideData("SolutionGradient");
dataCacheDouble &normals = cacheMapLeft.provideData("Normals");
dataCacheElement &cacheElementLeft = cacheMapLeft.getElement();
+ gradientSolutionLeft.set(fullMatrix<double>(group.getNbIntegrationPoints()*3,nbFields));
// required data
+ // inviscid
dataCacheDouble *boundaryTerm = boundaryCondition->newBoundaryTerm(cacheMapLeft);
- fullMatrix<double> normalFluxQP;
+ // viscous
+ dataCacheDouble *diffusiveFluxLeft = claw.newDiffusiveFlux(cacheMapLeft);
+ dataCacheDouble *maximumDiffusivityLeft = claw.newMaximumDiffusivity(cacheMapLeft);
+ dataCacheDouble *neumann = boundaryCondition->newDiffusiveNeumannBC(cacheMapLeft);
+ dataCacheDouble *dirichlet = boundaryCondition->newDiffusiveDirichletBC(cacheMapLeft);
+
+ fullMatrix<double> normalFluxQP,dPsiLeftDx,dofsLeft;
+
+ int p = group.getGroupLeft().getOrder();
+ int dim = group.getGroupLeft().getElement(0)->getDim();
for (int iFace=0 ; iFace<group.getNbElements() ;++iFace) {
normalFluxQP.setAsProxy(NormalFluxQP, iFace*nbFields, nbFields);
// ----- 2.3.1 --- provide the data to the cacheMap
solutionQPLeft.setAsProxy(solutionQP, iFace*nbFields, nbFields);
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);
+ dofsLeft.setAsProxy(solutionLeft, nbFields*group.getElementLeftId(iFace), nbFields);
+
uvw.setAsProxy(group.getIntegrationOnElementLeft(iFace));
cacheElementLeft.set(group.getElementLeft(iFace));
- // ----- 2.3.2 --- compute fluxes
+ // compute the gradient of the solution
+ if(gradientSolutionLeft.somethingDependOnMe()){
+ dPsiLeftDx.mult(dofsLeft, gradientSolutionLeft.set());
+ }
+
+ // ----- 2.3.2 --- compute inviscid contribution
for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
const double detJ = group.getDetJ (iFace, iPt);
for (int k=0;k<nbFields;k++)
normalFluxQP(iPt,k) = (*boundaryTerm)(iPt,k)*detJ;
}
+
+ if (diffusiveFluxLeft) {
+ for (int iPt =0; iPt< group.getNbIntegrationPoints(); iPt++) {
+ const double detJ = group.getDetJ (iFace, iPt);
+ //just for the lisibility :
+ 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;
+ 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);
delete boundaryTerm;
+ if (dirichlet) delete dirichlet;
+ if (neumann) delete neumann;
}
/*
@@ -520,6 +563,20 @@ 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
+// Then, groups of elements are created
+// -) 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 dgAlgorithm::buildGroups(GModel *model, int dim, int order,
std::vector<dgGroupOfElements*> &eGroups,
std::vector<dgGroupOfFaces*> &fGroups,
@@ -558,10 +615,15 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
allElements[entity->getMeshElement(iel)->getType()].push_back(entity->getMeshElement(iel));
}
}
+
+
// do a group of elements for every element type that is present in the mesh
Msg::Info("%d groups of elements",allElements.size());
- for (std::map<int, std::vector<MElement *> >::iterator it = allElements.begin(); it !=allElements.end() ; ++it)
+ for (std::map<int, std::vector<MElement *> >::iterator it = allElements.begin(); it !=allElements.end() ; ++it){
+
eGroups.push_back(new dgGroupOfElements(it->second,order));
+ }
+
for (int i=0;i<eGroups.size();i++){
// create a group of faces on the boundary for every group ef elements
@@ -612,40 +674,41 @@ void dgAlgorithm::residual( const dgConservationLaw &claw,
//volume term
for(size_t i=0;i<eGroups.size() ; i++) {
residu[i]->scale(0);
+
residualVolume(claw,*eGroups[i],*solution[i],*residu[i]);
}
// residu[0]->print("Volume");
//interface term
for(size_t i=0;i<fGroups.size() ; i++) {
- dgGroupOfFaces &faces = *fGroups[i];
- int iGroupLeft = -1, iGroupRight = -1;
- for(size_t j=0;j<eGroups.size() ; j++) {
- if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
- if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
- }
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
- faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
- residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
+ dgGroupOfFaces &faces = *fGroups[i];
+ int iGroupLeft = -1, iGroupRight = -1;
+ for(size_t j=0;j<eGroups.size() ; j++) {
+ if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+ if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
}
- // residu[0]->print("Interfaces");
- //boundaries
- for(size_t i=0;i<bGroups.size() ; i++) {
- dgGroupOfFaces &faces = *bGroups[i];
- int iGroupLeft = -1, iGroupRight = -1;
- for(size_t j=0;j<eGroups.size() ; j++) {
- if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
- if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
- }
- fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
- fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
- faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
- residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
- faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*2*nbFields);
+ faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
+ residualInterface(claw,faces,solInterface,*solution[iGroupLeft], *solution[iGroupRight],residuInterface);
+ faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
+ }
+ // residu[0]->print("Interfaces");
+ //boundaries
+ for(size_t i=0;i<bGroups.size() ; i++) {
+ dgGroupOfFaces &faces = *bGroups[i];
+ int iGroupLeft = -1, iGroupRight = -1;
+ for(size_t j=0;j<eGroups.size() ; j++) {
+ if (eGroups[j] == &faces.getGroupLeft())iGroupLeft = j;
+ if (eGroups[j] == &faces.getGroupRight())iGroupRight= j;
}
+ fullMatrix<double> solInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
+ fullMatrix<double> residuInterface(faces.getNbNodes(),faces.getNbElements()*nbFields);
+ faces.mapToInterface(nbFields, *solution[iGroupLeft], *solution[iGroupRight], solInterface);
+ residualBoundary(claw,faces,solInterface,*solution[iGroupLeft],residuInterface);
+ faces.mapFromInterface(nbFields, residuInterface, *residu[iGroupLeft], *residu[iGroupRight]);
+ }
- // residu[0]->print("Boundaries");
+ // residu[0]->print("Boundaries");
}
void dgAlgorithm::computeElementaryTimeSteps ( //dofManager &dof, // the DOF manager (maybe useless here)
diff --git a/Solver/dgConservationLaw.cpp b/Solver/dgConservationLaw.cpp
index 937866f44fb3fbad503bf9444137cc3ea279da0f..853b6756bbff9af99428e0e3c16a1f8599298010 100644
--- a/Solver/dgConservationLaw.cpp
+++ b/Solver/dgConservationLaw.cpp
@@ -1,23 +1,22 @@
#include "dgConservationLaw.h"
#include "function.h"
class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
- dgConservationLaw &_claw;
std::string _outsideValueFunctionName;
class term : public dataCacheDouble {
dataCacheMap cacheMapRight; // new cacheMap to pass to the Riemann solver
dataCacheDouble &solutionRight;
dataCacheDouble &outsideValue;
dataCacheDouble *riemannSolver;
- dgConservationLaw &_claw;
+ dgConservationLaw *_claw;
public:
- term(dgConservationLaw &claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
- dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw.nbFields()),
+ term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft,const std::string outsideValueFunctionName):
+ dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()),
cacheMapRight(cacheMapLeft.getNbEvaluationPoints()),
solutionRight(cacheMapRight.provideData("Solution")),
outsideValue(cacheMapLeft.get(outsideValueFunctionName,this)),
_claw(claw)
{
- riemannSolver=_claw.newRiemannSolver(cacheMapLeft,cacheMapRight);
+ riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapRight);
riemannSolver->addMeAsDependencyOf(this);
}
@@ -31,7 +30,7 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
}
};
public:
- dgBoundaryConditionOutsideValue(dgConservationLaw &claw,const std::string outsideValueFunctionName): _claw(claw),
+ dgBoundaryConditionOutsideValue(dgConservationLaw *claw,const std::string outsideValueFunctionName): dgBoundaryCondition(claw),
_outsideValueFunctionName(outsideValueFunctionName)
{ }
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
@@ -40,18 +39,17 @@ class dgBoundaryConditionOutsideValue : public dgBoundaryCondition {
};
class dgBoundarySymmetry : public dgBoundaryCondition {
- dgConservationLaw &_claw;
class term : public dataCacheDouble {
dataCacheDouble *riemannSolver;
- dgConservationLaw &_claw;
- public:
- term(dgConservationLaw &claw, dataCacheMap &cacheMapLeft):
- dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw.nbFields()), _claw(claw)
+ dgConservationLaw *_claw;
+ public:
+ term(dgConservationLaw *claw, dataCacheMap &cacheMapLeft):
+ dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()), _claw(claw)
{
- riemannSolver=_claw.newRiemannSolver(cacheMapLeft,cacheMapLeft);
+ riemannSolver=_claw->newRiemannSolver(cacheMapLeft,cacheMapLeft);
riemannSolver->addMeAsDependencyOf(this);
}
-
+
void _eval() {
if(riemannSolver){
for(int i=0;i<_value.size1(); i++)
@@ -60,39 +58,93 @@ class dgBoundarySymmetry : public dgBoundaryCondition {
}
}
};
- public:
- dgBoundarySymmetry(dgConservationLaw &claw): _claw(claw) {}
+public:
+ dgBoundarySymmetry(dgConservationLaw *claw): dgBoundaryCondition(claw) {}
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
return new term(_claw,cacheMapLeft);
}
};
class dgBoundaryCondition0Flux : public dgBoundaryCondition {
- dgConservationLaw &_claw;
class term : public dataCacheDouble {
- public:
- term(dgConservationLaw &claw,dataCacheMap &cacheMapLeft):
- dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw.nbFields()) {}
+ public:
+ term(dgConservationLaw *claw,dataCacheMap &cacheMapLeft):
+ dataCacheDouble(cacheMapLeft.getNbEvaluationPoints(),claw->nbFields()) {}
void _eval() {
}
};
- public:
- dgBoundaryCondition0Flux(dgConservationLaw &claw): _claw(claw) {}
+public:
+ dgBoundaryCondition0Flux(dgConservationLaw *claw): dgBoundaryCondition(claw) {}
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
return new term(_claw,cacheMapLeft);
}
};
dgBoundaryCondition *dgConservationLaw::newSymmetryBoundary() {
- return new dgBoundarySymmetry(*this);
+ return new dgBoundarySymmetry(this);
}
dgBoundaryCondition *dgConservationLaw::newOutsideValueBoundary(const std::string outsideValueFunctionName) {
- return new dgBoundaryConditionOutsideValue(*this,outsideValueFunctionName);
+ return new dgBoundaryConditionOutsideValue(this,outsideValueFunctionName);
}
dgBoundaryCondition *dgConservationLaw::new0FluxBoundary() {
- return new dgBoundaryCondition0Flux(*this);
+ return new dgBoundaryCondition0Flux(this);
}
+
+class dgBoundaryCondition::dirichlet_ : public dataCacheDouble {
+ dataCacheDouble /
+public:
+ dirichlet_(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this)){}
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,sol().size2());
+ for(int i=0; i< nQP; i++)
+ for (int k=0;k<sol().size2();k++)
+ _value(i,k) = sol(i,k);
+ }
+};
+
+class dgBoundaryCondition::neumann_ : public dataCacheDouble {
+ dgConservationLaw *_claw;
+ dataCacheDouble &sol,&normals;
+ dataCacheDouble *diffusiveFlux;
+public:
+ neumann_(dataCacheMap &cacheMap, dgConservationLaw *claw):
+ _claw (claw),
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)){
+ diffusiveFlux=_claw->newDiffusiveFlux(cacheMap);
+ if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
+ }
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,sol().size2());
+
+ const fullMatrix<double> &dfl = (*diffusiveFlux)();
+
+ for(int i=0; i< nQP; i++) {
+ for (int k=0;k<sol().size2();k++) {
+ _value(i,k) =
+ dfl(i,k+sol().size2()*0) *normals(0,i) +
+ dfl(i,k+sol().size2()*1) *normals(1,i) +
+ dfl(i,k+sol().size2()*2) *normals(2,i);
+ }
+ }
+ }
+ ~neumann_ () {if (diffusiveFlux)delete diffusiveFlux;}
+};
+
+dataCacheDouble *dgBoundaryCondition::newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const {
+ return new dirichlet_(cacheMapLeft);
+}
+dataCacheDouble *dgBoundaryCondition::newDiffusiveNeumannBC(dataCacheMap &cacheMapLeft) const {
+ return new dgBoundaryCondition::neumann_(cacheMapLeft, _claw);
+}
+
+
#include "Bindings.h"
void dgConservationLaw::registerBindings(binding *b){
diff --git a/Solver/dgConservationLaw.h b/Solver/dgConservationLaw.h
index 2bc50dc249558680184dc570257d51bdd5b5cd22..b93a2777846e6b15074650471c6f7b1295816e6d 100644
--- a/Solver/dgConservationLaw.h
+++ b/Solver/dgConservationLaw.h
@@ -14,9 +14,16 @@ class dgConservationLaw;
class binding;
class dgBoundaryCondition {
+ class neumann_;
+ class dirichlet_;
+ protected:
+ dgConservationLaw *_claw;
public:
+ dgBoundaryCondition (dgConservationLaw *claw) : _claw(claw){}
virtual ~dgBoundaryCondition () {}
virtual dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const = 0;
+ virtual dataCacheDouble *newDiffusiveNeumannBC(dataCacheMap &cacheMapLeft) const;
+ virtual dataCacheDouble *newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const;
static void registerBindings(binding *b);
};
diff --git a/Solver/dgConservationLawPerfectGas.cpp b/Solver/dgConservationLawPerfectGas.cpp
index 8e221b7e4ecaf24be1064d0f79ed88e937b8858f..9160134f075f33120a58ad0fd4ee22540123b53a 100644
--- a/Solver/dgConservationLawPerfectGas.cpp
+++ b/Solver/dgConservationLawPerfectGas.cpp
@@ -1,6 +1,10 @@
#include "dgConservationLawPerfectGas.h"
#include "function.h"
+/*
+\rho \left(\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \left( \tfrac13 \mu + \mu^v) \nabla (\nabla \cdot \mathbf{v} \right) + \rho \mathbf{f}
+*/
+
const double GAMMA = 1.4;
static inline void _ROE2D (const double &_GAMMA,
const double &nx,
@@ -169,6 +173,89 @@ class dgPerfectGasLaw2d::advection : public dataCacheDouble {
}
};
+class dgPerfectGasLaw2d::diffusion : public dataCacheDouble {
+ dataCacheDouble &sol,&grad,&mu,κ
+ public:
+ diffusion(dataCacheMap &cacheMap,
+ const std::string &muFunctionName,
+ const std::string &kappaFunctionName ):
+ sol(cacheMap.get("Solution",this)),
+ grad(cacheMap.get("SolutionGradient",this)),
+ mu(cacheMap.get(muFunctionName,this)),
+ kappa(cacheMap.get(kappaFunctionName,this))
+ {};
+ void _eval () {
+ const int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value=fullMatrix<double>(nQP,12);
+
+ for (size_t k = 0 ; k < nQP; k++ ){
+ const double muk = mu(k,0);
+ const double kappak = kappa(k,0);
+ // printf("%g %g grad rho (%g %g)\n",muk,kappak,grad(k*3+0,0),grad(k*3+1,0));
+ // find gradients of primitive variables
+ const double overRho = 1./sol(k,0);
+ const double u = sol(k,1) * overRho;
+ const double v = sol(k,2) * overRho;
+ // grad v
+ const double dudx = (grad(k*3+0,1) - grad(k*3+0,0) * u) * overRho;
+ const double dudy = (grad(k*3+1,1) - grad(k*3+1,0) * u) * overRho;
+ const double dvdx = (grad(k*3+0,2) - grad(k*3+0,0) * v) * overRho;
+ const double dvdy = (grad(k*3+1,2) - grad(k*3+1,0) * v) * overRho;
+ // shear stress tensor - Newtonian fluid without bulk viscosity
+
+ const double tauxx = muk * (dudx - dvdy);
+ const double tauxy = muk * (dudy + dvdx);
+ const double tauyy = - tauxx;
+
+ _value(k,1) = -tauxx;
+ _value(k,2) = -tauxy;
+
+ _value(k,1+4) = -tauxy;
+ _value(k,2+4) = -tauyy;
+ // heat flux
+
+ const double kOverRhoCv = kappak * overRho;
+ const double Ek = u * u + v * v;
+ const double eMinus= sol(k,3) * overRho - Ek;
+
+ // k grad T
+
+ _value(k,3) = -kOverRhoCv*(grad(k*3+0,3) - eMinus*grad(k*3+0,0) - (grad(k*3+0,1)*u + grad(k*3+0,2)*v));
+ _value(k,3+4) = -kOverRhoCv*(grad(k*3+1,3) - eMinus*grad(k*3+1,0) - (grad(k*3+1,1)*u + grad(k*3+1,2)*v));
+
+ // v . tau - momentum dissipation
+
+ _value(k,3) -= (u * tauxx + v * tauxy);
+ _value(k,3+4) -= (u * tauxy + v * tauyy);
+
+ // printf("%g %g %g %g %g %g %g %g \n",_value(k,0),_value(k,1),_value(k,2),_value(k,3),_value(k,4),_value(k,5),_value(k,6),_value(k,7));
+
+ }
+ }
+};
+
+class dgPerfectGasLaw2d::source : public dataCacheDouble {
+ dataCacheDouble &sol,&s;
+public:
+ source(dataCacheMap &cacheMap, const std::string &sourceFunctionName):
+ sol(cacheMap.get("Solution",this)),
+ s(cacheMap.get(sourceFunctionName,this))
+ {};
+ void _eval () {
+ const int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value=fullMatrix<double>(nQP,4);
+ for (size_t k = 0 ; k < nQP; k++ ){
+ _value(k,0) = sol(k,0)*s(0,0);
+ _value(k,1) = sol(k,0)*s(0,1);
+ _value(k,2) = sol(k,0)*s(0,2);
+ _value(k,3) = sol(k,0)*s(0,3);
+ }
+ }
+};
+
+
class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
dataCacheDouble &normals, &solL, &solR;
public:
@@ -180,7 +267,7 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
void _eval () {
int nQP = solL().size1();
if(_value.size1() != nQP)
- _value = fullMatrix<double>(nQP,8);
+ _value = fullMatrix<double>(nQP,12);
for(int i=0; i< nQP; i++) {
const double solLeft [4] = {solL(i,0),solL(i,1),solL(i,2),solL(i,3)};
@@ -201,6 +288,55 @@ class dgPerfectGasLaw2d::riemann : public dataCacheDouble {
}
}
};
+
+class dgPerfectGasLaw2d::maxConvectiveSpeed : public dataCacheDouble {
+ dataCacheDouble /
+ public:
+ maxConvectiveSpeed(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this))
+ {
+ };
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value=fullMatrix<double>(nQP,1);
+ for(int i=0; i< nQP; i++) {
+ const double rhov2 = (sol(i,1)*sol(i,1) + sol(i,2)*sol(i,2))/sol(i,0);
+ const double p = (GAMMA-1.0)*(sol(i,3) - 0.5 * rhov2);
+ // printf("p = %g %g %g %g %g\n",p,_value(i,0),_value(i,1),_value(i,2),_value(i,3));
+ const double c = sqrt(GAMMA*p/sol(i,0));
+ _value(i,0) = c+sqrt(rhov2/sol(i,0));
+ }
+ }
+};
+
+class dgPerfectGasLaw2d::maxDiffusivity : public dataCacheDouble {
+ dataCacheDouble &sol,&mu,κ
+ public:
+ maxDiffusivity(dataCacheMap &cacheMap, const std::string &muFunctionName, const std::string &kappaFunctionName ):
+ sol(cacheMap.get("Solution",this)),
+ mu(cacheMap.get(muFunctionName,this)),
+ kappa(cacheMap.get(kappaFunctionName,this))
+ {
+ };
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value=fullMatrix<double>(nQP,1);
+ for(int i=0; i< nQP; i++) {
+ _value(i,0) = std::max(mu(i,0),kappa(i,0));
+ }
+ }
+};
+
+
+dataCacheDouble *dgPerfectGasLaw2d::newMaximumDiffusivity( dataCacheMap &cacheMap) const {
+ return _muFunctionName.empty() ? NULL : new maxDiffusivity(cacheMap,_muFunctionName,_kappaFunctionName);
+}
+
+dataCacheDouble *dgPerfectGasLaw2d::newMaxConvectiveSpeed( dataCacheMap &cacheMap) const {
+ return new maxConvectiveSpeed(cacheMap);
+}
dataCacheDouble *dgPerfectGasLaw2d::newConvectiveFlux( dataCacheMap &cacheMap) const {
return new advection(cacheMap);
}
@@ -208,24 +344,38 @@ dataCacheDouble *dgPerfectGasLaw2d::newRiemannSolver( dataCacheMap &cacheMapLeft
return new riemann(cacheMapLeft, cacheMapRight);
}
dataCacheDouble *dgPerfectGasLaw2d::newDiffusiveFlux( dataCacheMap &cacheMap) const {
- return 0;
+ if (_muFunctionName.empty() || _kappaFunctionName.empty())
+ return 0;
+ else
+ return new diffusion(cacheMap,_muFunctionName,_kappaFunctionName);
}
dataCacheDouble *dgPerfectGasLaw2d::newSourceTerm (dataCacheMap &cacheMap) const {
- return 0;
+ if (_sourceFunctionName.empty())
+ return 0;
+ else
+ return new source(cacheMap,_sourceFunctionName);
}
+
dgPerfectGasLaw2d::dgPerfectGasLaw2d()
{
_nbf = 4; // \rho \rho u \rho v \rho e
}
-
+//-------------------------------------------------------------------------------
+// A class for slip and non slip walls
+// could easily add moving walls ...
+//-------------------------------------------------------------------------------
class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
+ int _type;
class term : public dataCacheDouble {
dataCacheDouble &sol,&normals;
- public:
+ public:
+//-------------------------------------------------------------------------------
+// NON VISCOUS BOUNDARY FLUX
+//-------------------------------------------------------------------------------
term(dataCacheMap &cacheMap):
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)){}
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)){}
void _eval () {
int nQP = sol().size1();
if(_value.size1() != nQP)
@@ -240,8 +390,13 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
sol(i,1) - 2 * vn * nx,
sol(i,2) - 2 * vn * ny,
sol(i,3)};
- // double FLUX[4] ;
- // _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX);
+ double FLUX[4] ;
+ _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX);
+ _value(i,0) = FLUX[0];
+ _value(i,1) = FLUX[1];
+ _value(i,2) = FLUX[2];
+ _value(i,3) = FLUX[3];
+ /*
const double q11 = sol(i,1)*sol(i,1)/sol(i,0);
const double q22 = sol(i,2)*sol(i,2)/sol(i,0);
const double p = (GAMMA-1)*sol(i,3) - 0.5*(GAMMA-1)*(q11+q22);
@@ -249,72 +404,153 @@ class dgBoundaryConditionPerfectGasLaw2dWall : public dgBoundaryCondition {
_value(i,1) = -p*nx;//FLUX[1];
_value(i,2) = -p*ny;//FLUX[2];
_value(i,3) = 0.0;//FLUX[3];
+ */
}
}
};
- public:
- dgBoundaryConditionPerfectGasLaw2dWall(){}
- dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
- return new term(cacheMapLeft);
- }
-};
-dgBoundaryCondition *dgPerfectGasLaw2d::newWallBoundary() const {
- return new dgBoundaryConditionPerfectGasLaw2dWall();
-}
-#if 0 // JF : I commented this out as I think this equivalent to the generic OutsideValue condition
-can you confirm and remove it ?
-class dgBoundaryConditionPerfectGasLaw2dFreeStream : public dgBoundaryCondition {
- class term : public dataCacheDouble {
- dataCacheDouble &sol,&normals,&freeStream;
+//-------------------------------------------------------------------------------
+// Non Slip Dirichlet BC --
+// What to put on the right side of the IP (no velocity, same density)
+// Assume Adiabatic Wall --> no heat transfer -> \partial_n T = 0
+// So, we use the same temperature as inside
+//-------------------------------------------------------------------------------
+ class dirichletNonSlip : public dataCacheDouble {
+ dataCacheDouble /
public:
- term(dataCacheMap &cacheMap, const std::string & freeStreamName):
- sol(cacheMap.get("Solution",this)),
- normals(cacheMap.get("Normals",this)),
- freeStream(cacheMap.get(freeStreamName,this)){}
+ dirichletNonSlip(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this)){}
void _eval () {
int nQP = sol().size1();
if(_value.size1() != nQP)
_value = fullMatrix<double>(nQP,4);
for(int i=0; i< nQP; i++) {
- const double nx = normals(0,i);
- const double ny = normals(1,i);
- const double solLeft [4] = {sol(i,0),sol(i,1),sol(i,2),sol(i,3)};
- const double solRight[4] = {freeStream(i,0),
- freeStream(i,1),
- freeStream(i,2),
- freeStream(i,3)};
- double FLUX[4] ;
- _ROE2D (GAMMA,nx,ny,solLeft,solRight,FLUX);
- /*
- printf("SOLL %g %g %g %g\n",solLeft[0],solLeft[1],solLeft[2], solLeft[3]);
- printf("SOLR %g %g %g %g\n",solRight[0],solRight[1],solRight[2], solRight[3]);
- printf("N %g %g\n",nx,ny);
- printf("FLUX %g %g %g %g\n",FLUX[0],FLUX[1],FLUX[2],FLUX[3]);
- */
- _value(i,0) = FLUX[0];
- _value(i,1) = FLUX[1];
- _value(i,2) = FLUX[2];
- _value(i,3) = FLUX[3];
+ _value(i,0) = sol(i,0);
+ _value(i,1) = 0.0;
+ _value(i,2) = 0.0;
+ _value(i,3) = sol(i,3) - 0.5 * (sol(i,1) * sol(i,1) + sol(i,2) * sol(i,2))/sol(i,0);
}
}
};
+
+//-------------------------------------------------------------------------------
+// Non Slip Neumann BC --
+// Compute normal diffusive fluxes at the boundary
+// Assume Adiabatic Wall --> no heat transfer -> no flux for component 3 (in 2D)
+//-------------------------------------------------------------------------------
+
+ class neumannNonSlip : public dataCacheDouble {
+ dgConservationLaw *_claw;
+ dataCacheDouble &sol,&normals;
+ dataCacheDouble *diffusiveFlux;
public:
- std::string _freeStreamName;
- dgBoundaryConditionPerfectGasLaw2dFreeStream(std::string & freeStreamName)
- : _freeStreamName(freeStreamName){}
- dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
- return new term(cacheMapLeft,_freeStreamName);
- }
+ neumannNonSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
+ _claw (claw),
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)){
+ diffusiveFlux=_claw->newDiffusiveFlux(cacheMap);
+ if (diffusiveFlux)diffusiveFlux->addMeAsDependencyOf(this);
+ }
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,4);
+
+ const fullMatrix<double> &dfl = (*diffusiveFlux)();
+
+ for(int i=0; i< nQP; i++) {
+ for (int k=0;k<3;k++) {
+ _value(i,k) =
+ dfl(i,k+4*0) *normals(0,i) +
+ dfl(i,k+4*1) *normals(1,i) +
+ dfl(i,k+4*2) *normals(2,i);
+ }
+ _value(i,3) = 0.0;
+ }
+ }
+ ~neumannNonSlip () {if (diffusiveFlux)delete diffusiveFlux;}
+ };
+
+//-------------------------------------------------------------------------------
+// Slip Wall Dirichlet BC --
+// Assume zero derivatives of all variables --> put the same as inside
+//-------------------------------------------------------------------------------
+
+ class dirichletSlip : public dataCacheDouble {
+ dataCacheDouble /
+ public:
+ dirichletSlip(dataCacheMap &cacheMap):
+ sol(cacheMap.get("Solution",this)){}
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,sol().size2());
+
+ for(int i=0; i< nQP; i++) {
+ for(int k=0; k< sol().size2(); k++) {
+ _value(i,k) = sol(i,k);
+ }
+ }
+ }
+ };
+
+//-------------------------------------------------------------------------------
+// Slip Wall or Symmetry Neumann BC -- assume NO FLUXES AT ALL
+//-------------------------------------------------------------------------------
+
+ class neumannSlip : public dataCacheDouble {
+ dgConservationLaw *_claw;
+ dataCacheDouble &sol,&normals;
+ public:
+ neumannSlip(dataCacheMap &cacheMap, dgConservationLaw *claw):
+ _claw (claw),
+ sol(cacheMap.get("Solution",this)),
+ normals(cacheMap.get("Normals",this)){
+ }
+ void _eval () {
+ int nQP = sol().size1();
+ if(_value.size1() != nQP)
+ _value = fullMatrix<double>(nQP,4);
+ _value.setAll(0.0);
+ }
+ };
+
+ public:
+ dgBoundaryConditionPerfectGasLaw2dWall(int type, dgPerfectGasLaw2d *claw):_type(type), dgBoundaryCondition(claw){}
+ dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
+ return new term(cacheMapLeft);
+ }
+ dataCacheDouble *newDiffusiveDirichletBC(dataCacheMap &cacheMapLeft) const {
+ switch(_type){
+ case 1 : return new dirichletSlip(cacheMapLeft);
+ case 2 : return new dirichletNonSlip(cacheMapLeft);
+ }
+ }
+ dataCacheDouble *newDiffusiveNeumannBC(dataCacheMap &cacheMapLeft) const {
+ switch(_type){
+ case 1 : return new neumannSlip(cacheMapLeft,_claw);
+ case 2 : return new neumannNonSlip(cacheMapLeft,_claw);
+ }
+ }
};
-#endif
+
+dgBoundaryCondition *dgPerfectGasLaw2d::newWallBoundary() {
+ return new dgBoundaryConditionPerfectGasLaw2dWall(2, this);
+}
+dgBoundaryCondition *dgPerfectGasLaw2d::newSlipWallBoundary() {
+ return new dgBoundaryConditionPerfectGasLaw2dWall(1, this);
+}
#include "Bindings.h"
void dgPerfectGasLaw2dRegisterBindings (binding *b){
classBinding *cb = b->addClass<dgPerfectGasLaw2d>("dgPerfectGasLaw2d");
methodBinding *cm;
cb->addMethod("newWallBoundary",&dgPerfectGasLaw2d::newWallBoundary);
+ cb->addMethod("newNonSlipWallBoundary",&dgPerfectGasLaw2d::newWallBoundary);
+ cb->addMethod("newSlipWallBoundary",&dgPerfectGasLaw2d::newSlipWallBoundary);
+ cb->addMethod("setSource",&dgPerfectGasLaw2d::setSource);
+ cb->addMethod("setViscosityAndThermalConductivity",&dgPerfectGasLaw2d::setViscosityAndThermalConductivity);
cb->setConstructor<dgPerfectGasLaw2d>();
cb->setParentClass<dgConservationLaw>();
}
diff --git a/Solver/dgConservationLawPerfectGas.h b/Solver/dgConservationLawPerfectGas.h
index c53ff8c057d001a3b642873c3d5ad7f8bd402e2d..75e43d9a2c716868b1017da666c118780b462cbc 100644
--- a/Solver/dgConservationLawPerfectGas.h
+++ b/Solver/dgConservationLawPerfectGas.h
@@ -2,16 +2,37 @@
#define _DG_CONSERVATION_LAW_PERFECT_GAS_H_
#include "dgConservationLaw.h"
class dgPerfectGasLaw2d : public dgConservationLaw {
- // perfect gas law, GAMMA is the only parameter
class advection;
+ class diffusion;
class riemann;
+ class source;
+ class maxConvectiveSpeed;
+ class maxDiffusivity;
+ // the name of the functions for
+ // viscosity (_muFunctionName)
+ // thermal conductivity (_kappaFunctionName)
+ std::string _kappaFunctionName,_muFunctionName;
+ std::string _sourceFunctionName;
+
public:
+ dataCacheDouble *newMaxConvectiveSpeed( dataCacheMap &cacheMap) const;
+ dataCacheDouble *newMaximumDiffusivity( dataCacheMap &cacheMap) const;
dataCacheDouble *newConvectiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newRiemannSolver( dataCacheMap &cacheMapLeft, dataCacheMap &cacheMapRight) const;
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const;
- dgPerfectGasLaw2d();
- dgBoundaryCondition *newWallBoundary()const ;
+ dgPerfectGasLaw2d();
+ dgBoundaryCondition *newWallBoundary() ;
+ dgBoundaryCondition *newSlipWallBoundary() ;
+ // sets some coefficients
+ void setViscosityAndThermalConductivity (std::string a, std::string b){
+ _muFunctionName = a;
+ _kappaFunctionName = b;
+ }
+ // sets some coefficients
+ void setSource(std::string a){
+ _sourceFunctionName = a;
+ }
};
class binding;
void dgPerfectGasLaw2dRegisterBindings(binding *b);
diff --git a/Solver/dgConservationLawShallowWater2d.cpp b/Solver/dgConservationLawShallowWater2d.cpp
index 5a1e4b68d132571db7aec41c6a268eacaf3ed126..b0fe3019c2208778b8f8d661ccc62765067d136c 100644
--- a/Solver/dgConservationLawShallowWater2d.cpp
+++ b/Solver/dgConservationLawShallowWater2d.cpp
@@ -117,7 +117,7 @@ class dgConservationLawShallowWater2d::boundaryWall : public dgBoundaryCondition
}
};
public:
- boundaryWall(){}
+ boundaryWall(dgConservationLaw *claw) : dgBoundaryCondition(claw){}
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
return new term(cacheMapLeft);
}
@@ -137,7 +137,7 @@ dataCacheDouble *dgConservationLawShallowWater2d::newSourceTerm (dataCacheMap &c
}
dgBoundaryCondition *dgConservationLawShallowWater2d::newBoundaryWall(){
- return new boundaryWall();
+ return new boundaryWall(this);
}
#include "Bindings.h"
diff --git a/Solver/dgConservationLawWaveEquation.cpp b/Solver/dgConservationLawWaveEquation.cpp
index 893da29c76f7ebd86480fc9d12508d5b19c48183..e7c27a94b790c8e6032bb37aa822477a6fd7335b 100644
--- a/Solver/dgConservationLawWaveEquation.cpp
+++ b/Solver/dgConservationLawWaveEquation.cpp
@@ -126,13 +126,13 @@ class dgBoundaryConditionWaveEquationWall : public dgBoundaryCondition {
}
};
public:
- dgBoundaryConditionWaveEquationWall(int DIM):_DIM(DIM){}
+ dgBoundaryConditionWaveEquationWall(int DIM, dgConservationLaw *claw):_DIM(DIM),dgBoundaryCondition(claw){}
dataCacheDouble *newBoundaryTerm(dataCacheMap &cacheMapLeft) const {
return new term(cacheMapLeft,_DIM);
}
};
-dgBoundaryCondition *dgConservationLawWaveEquation::newBoundaryWall()const{
- return new dgBoundaryConditionWaveEquationWall(_DIM);
+dgBoundaryCondition *dgConservationLawWaveEquation::newBoundaryWall(){
+ return new dgBoundaryConditionWaveEquationWall(_DIM, this);
}
#include "Bindings.h"
diff --git a/Solver/dgConservationLawWaveEquation.h b/Solver/dgConservationLawWaveEquation.h
index ee3498acdd00ded72aecc1712a91e7c0b7c2ab38..0202905d09a34ffa14a68b8d23b316be2cc33e04 100644
--- a/Solver/dgConservationLawWaveEquation.h
+++ b/Solver/dgConservationLawWaveEquation.h
@@ -14,7 +14,7 @@ class dgConservationLawWaveEquation : public dgConservationLaw {
dataCacheDouble *newDiffusiveFlux( dataCacheMap &cacheMap) const ;
dataCacheDouble *newSourceTerm (dataCacheMap &cacheMap) const ;
dataCacheDouble *newMaxConvectiveSpeed (dataCacheMap &cacheMap) const;
- dgBoundaryCondition *newBoundaryWall () const;
+ dgBoundaryCondition *newBoundaryWall () ;
dgConservationLawWaveEquation(int DIM);
static const char *className;
static const char *parentClassName;
diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp
index 211c9c684428eba7285a85032f97e1d3ed199a56..f8350962bb529f5a9f6d8eb34f32ce62ba74f36f 100644
--- a/Solver/dgSystemOfEquations.cpp
+++ b/Solver/dgSystemOfEquations.cpp
@@ -69,6 +69,8 @@ void dgSystemOfEquations::registerBindings(binding *b){
cb->addMethod("computeInvSpectralRadius",&dgSystemOfEquations::computeInvSpectralRadius);
cb->addMethod("RK44_limiter",&dgSystemOfEquations::RK44_limiter);
cb->addMethod("multirateRK43",&dgSystemOfEquations::multirateRK43);
+ cb->addMethod("saveSolution",&dgSystemOfEquations::saveSolution);
+ cb->addMethod("loadSolution",&dgSystemOfEquations::loadSolution);
}
// do a L2 projection
@@ -143,6 +145,18 @@ dgSystemOfEquations::~dgSystemOfEquations(){
}
}
+void dgSystemOfEquations::saveSolution (const std::string &name) const{
+ FILE *f = fopen (name.c_str(),"wb");
+ _solution->_data->binarySave(f);
+ fclose(f);
+}
+
+void dgSystemOfEquations::loadSolution (const std::string &name){
+ FILE *f = fopen (name.c_str(),"rb");
+ _solution->_data->binaryLoad(f);
+ fclose(f);
+}
+
void dgSystemOfEquations::export_solution_as_is (const std::string &name){
// the elementnodedata::export does not work !!
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index edc36dddf71bafa1f5ede84d035e2e6c01800939..5799445d4d253de7cd3aa0ac7b15fa17e2a14282 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -71,6 +71,8 @@ public:
iElement * _claw->nbFields(),_claw->nbFields());
}
void export_solution_as_is (const std::string &fileName);
+ void saveSolution (const std::string &fileName) const;
+ void loadSolution (const std::string &fileName);
~dgSystemOfEquations();
private:
};
diff --git a/benchmarks/2d/square.geo b/benchmarks/2d/square.geo
index 88dbb82010bc3abe731410df352f0df65a5b9c95..4ef669b84647aa5c773b28f5cc0c6e2212ecfe6f 100644
--- a/benchmarks/2d/square.geo
+++ b/benchmarks/2d/square.geo
@@ -11,5 +11,5 @@ Line(4) = {1, 2};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(10) = {5};
-Compound Surface(-11)={10};
+//Compound Surface(-11)={10};