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

dg : correct a few small bugs in lua bindings and replace main files by lua input files

parent 6afba50e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 100 additions and 339 deletions
CMakeLists.txt
+ 0
− 8
View file @ b4ffa082
......@@ -1005,16 +1005,8 @@ message("")
mark_as_advanced(BISON FLEX GMP_LIB GMSH_EXTRA_VERSION HDF5_LIB MAKEINFO
MED_LIB OCC_INC SZ_LIB TAUCS_LIB LUA_LIB TEXI2PDF)
add_executable(gmshdg EXCLUDE_FROM_ALL Solver/dgMainTest.cpp ${GMSH_SRC})
add_executable(gmshlua EXCLUDE_FROM_ALL Solver/dgMainLua.cpp ${GMSH_SRC})
target_link_libraries(gmshdg ${LINK_LIBRARIES})
target_link_libraries(gmshlua ${LINK_LIBRARIES})
add_executable(gmshdgsw EXCLUDE_FROM_ALL Solver/dgMainShallowWater2d.cpp ${GMSH_SRC})
target_link_libraries(gmshdgsw ${LINK_LIBRARIES})
add_executable(gmshdgwave EXCLUDE_FROM_ALL Solver/dgMainWaveEquation.cpp ${GMSH_SRC})
target_link_libraries(gmshdgwave ${LINK_LIBRARIES})
add_executable(gmshdgsam EXCLUDE_FROM_ALL Solver/dgMainSam.cpp ${GMSH_SRC})
target_link_libraries(gmshdgsam ${LINK_LIBRARIES})
Solver/TESTCASES/Advection.lua 0 → 100644
+ 45
− 0
View file @ b4ffa082
model = GModel ()
model:load ('square.geo')
model:load ('square.msh')
dg = dgSystemOfEquations (model)
dg:setOrder(1)
-- conservation law
-- advection speed
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))
-- boundary condition
outside=fullMatrix(1,1)
outside:set(0,0,0.)
dg:addBoundaryCondition('Border','OutsideValues',createFunction.constant(outside))
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/Advection_00000')
-- main loop
for i=1,800 do
norm = dg:RK44(0.01)
if (i % 10 == 0) then
print('iter',i,norm)
dg:exportSolution(string.format("output/Advection-%05d", i))
end
end
Solver/TESTCASES/Stommel.lua 0 → 100644
+ 18
− 0
View file @ b4ffa082
model = GModel()
model:load ('stommel_square.msh')
dg = dgSystemOfEquations (model)
dg:setOrder (1)
dg:setConservationLaw('ShallowWater2d')
dg:addBoundaryCondition('Wall','Wall')
dg:setup()
dg:exportSolution('output/solution_00000')
for i=1,1000 do
norm = dg:RK44(50)
if ( i%100 ==0 ) then
print ('iter ', i, norm)
dg:exportSolution(string.format('output/solution-%05d',i))
end
end
Solver/TESTCASES/WavePulse.lua
+ 2
− 2
View file @ b4ffa082
......@@ -39,7 +39,7 @@ DG:L2Projection(initialCondition)
print'*** export ***'
DG:exportSolution('solution_0')
DG:exportSolution('output/solution_000')
print'*** solve ***'
......@@ -48,7 +48,7 @@ for i=1,1000 do
norm = DG:RK44(dt)
print('*** ITER ***',i,norm)
if (i % 10 == 0) then
AA = string.format("solution-%03d", i)
AA = string.format("output/solution-%03d", i)
DG:exportSolution(AA)
end
end
......
Solver/TESTCASES/stommel_square.geo 0 → 100644
+ 16
− 0
View file @ b4ffa082
l = 5e5;
Point(1) = {-l, -l, 0};
Point(2) = {l, -l, 0};
Point(3) = {l, l, 0};
Point(4) = {-l, l, 0};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {3, 4, 1, 2};
Field[1] = MathEval;
Field[1].F = "0.5e5";
Background Field = 1;
Plane Surface(6) = {5};
Physical Surface(0)={6};
Physical Line("Wall") = {3, 2, 1, 4};
Solver/dgMainShallowWater2d.cpp deleted 100644 → 0
+ 0
− 71
View file @ 6afba50e
#include <stdio.h>
#include <vector>
#include "GModel.h"
#include "dgGroupOfElements.h"
#include "dgAlgorithm.h"
#include "dgConservationLaw.h"
#include "Gmsh.h"
#include "function.h"
#include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v,int iField=1, int nbFields=1);
int main(int argc, char **argv) {
GmshMergeFile("input/stommel.msh");
//we probably need a class to group those three ones
std::vector<dgGroupOfElements*> elementGroups;
std::vector<dgGroupOfFaces*> faceGroups;
std::vector<dgGroupOfFaces*> boundaryGroups;
int order=1;
int dimension=2;
dgAlgorithm algo;
function::registerDefaultFunctions();
algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
//for now, we suppose there is only one group of elements
int nbNodes = elementGroups[0]->getNbNodes();
fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*3);
fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()*3);
dgConservationLaw *law = dgNewConservationLawShallowWater2d();
law->addBoundaryCondition("Wall",dgNewBoundaryConditionShallowWater2dWall());
for(int iT=0; iT<10000; iT++) {
if(iT%100==0){
printf("%i\n",iT);
char name[100];
sprintf(name,"output/eta_%05i.pos",iT/100);
print(name,*elementGroups[0],&sol(0,0),0,3);
sprintf(name,"output/u_%05i.pos",iT/100);
print(name,*elementGroups[0],&sol(0,0),1,3);
sprintf(name,"output/v_%05i.pos",iT/100);
print(name,*elementGroups[0],&sol(0,0),2,3);
}
algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,50,residu,sol);
}
delete law;
}
void print (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields) {
FILE *file = fopen(filename,"w");
fprintf(file,"View \"%s\" {\n", filename);
for(int iel=0;iel<els.getNbElements();iel++){
MElement *el = els.getElement(iel);
fprintf(file,"ST (");
for (int iv=0; iv<el->getNumVertices(); iv++) {
MVertex *vertex = el->getVertex(iv);
SPoint3 coord = vertex->point();
fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
}
fprintf(file,"{");
v+=(iField)*el->getNumVertices();
for (int iv=0; iv<el->getNumVertices(); iv++){
fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':',');
}
v+=(nbFields-1-iField)*el->getNumVertices();
fprintf(file,";\n");
}
fprintf(file,"};");
fclose(file);
}
Solver/dgMainTest.cpp deleted 100644 → 0
+ 0
− 111
View file @ 6afba50e
#include <stdio.h>
#include <vector>
#include "GModel.h"
#include "dgGroupOfElements.h"
#include "dgAlgorithm.h"
#include "dgConservationLaw.h"
#include "Gmsh.h"
#include "function.h"
#include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v);
class dgConservationLawInitialCondition : public dgConservationLaw {
class gaussian : public dataCacheDouble {
dataCacheDouble &xyz;
double _xc,_yc;
public:
gaussian(dataCacheMap &cacheMap,double xc, double yc):xyz(cacheMap.get("XYZ",this)),_xc(xc),_yc(yc){};
void _eval () {
if(_value.size1() != xyz().size1())
_value=fullMatrix<double>(xyz().size1(),1);
for(int i=0; i< _value.size1(); i++) {
_value(i,0)=exp(-(pow(xyz(i,0)-_xc,2)+pow(xyz(i,1)-_yc,2))*100);
}
}
};
public:
dgConservationLawInitialCondition() {
_nbf = 1;
}
dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
return new gaussian(cacheMap,0.2,0.3);
}
};
int main(int argc, char **argv){
GmshMergeFile("input/square1.msh");
//we probably need a class to group those three ones
std::vector<dgGroupOfElements*> elementGroups;
std::vector<dgGroupOfFaces*> faceGroups;
std::vector<dgGroupOfFaces*> boundaryGroups;
int order=1;
int dimension=2;
dgAlgorithm algo;
function::registerDefaultFunctions();
algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
//for now, we suppose there is only one group of elements
int nbNodes = elementGroups[0]->getNbNodes();
fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements());
fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements());
// initial condition
{
dgConservationLawInitialCondition initLaw;
algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
}
fullMatrix<double> advectionSpeed(3,1);
advectionSpeed(0,0)=0.15;
advectionSpeed(1,0)=0.05;
advectionSpeed(2,0)=0.;
function::add("advectionSpeed",function::newFunctionConstant(advectionSpeed));
dgConservationLaw *law = dgNewConservationLawAdvection("advectionSpeed");
fullMatrix<double> outsideValue(1,1);
outsideValue(0,0)=0;
function::add("outsideValue",function::newFunctionConstant(outsideValue));
law->addBoundaryCondition("Left",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue"));
law->addBoundaryCondition("Right",dgBoundaryCondition::newOutsideValueCondition(*law,"outsideValue"));
law->addBoundaryCondition("Top",dgBoundaryCondition::new0FluxCondition(*law));
law->addBoundaryCondition("Bottom",dgBoundaryCondition::new0FluxCondition(*law));
print("output/init.pos",*elementGroups[0],&sol(0,0));
for(int iT=0; iT<100; iT++) {
algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,0.01,residu,sol);
if(iT%10==0){
char name[100];
sprintf(name,"output/test_%05i.pos",iT/10);
printf("%i\n",iT);
print(name,*elementGroups[0],&sol(0,0));
}
}
delete law;
}
void print (const char *filename,const dgGroupOfElements &els, double *v) {
FILE *file = fopen(filename,"w");
fprintf(file,"View \"%s\" {\n", filename);
for(int iel=0;iel<els.getNbElements();iel++){
MElement *el = els.getElement(iel);
fprintf(file,"ST (");
for (int iv=0; iv<el->getNumVertices(); iv++) {
MVertex *vertex = el->getVertex(iv);
SPoint3 coord = vertex->point();
fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
}
fprintf(file,"{");
for (int iv=0; iv<el->getNumVertices(); iv++)
fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':',');
fprintf(file,";\n");
}
fprintf(file,"};");
fclose(file);
}
Solver/dgMainWaveEquation.cpp deleted 100644 → 0
+ 0
− 141
View file @ 6afba50e
#include <stdio.h>
#include <vector>
#include "GModel.h"
#include "dgGroupOfElements.h"
#include "dgAlgorithm.h"
#include "dgConservationLaw.h"
#include "Gmsh.h"
#include "function.h"
#include "functionDerivator.h"
#include "MElement.h"
void print (const char *filename,const dgGroupOfElements &els, double *v,int iField=1, int nbFields=1);
void print_vector (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields);
class dgConservationLawInitialCondition : public dgConservationLaw {
class gaussian : public dataCacheDouble {
dataCacheDouble &xyz;
double _xc,_yc;
public:
gaussian(dataCacheMap &cacheMap,double xc, double yc):xyz(cacheMap.get("XYZ",this)),_xc(xc),_yc(yc){};
void _eval () {
if(_value.size1() != xyz().size1())
_value=fullMatrix<double>(xyz().size1(),3);
for(int i=0; i< _value.size1(); i++) {
_value(i,0)=exp(-(pow(xyz(i,0)-_xc,2)+pow(xyz(i,1)-_yc,2))*100)+1;
}
}
};
public:
dgConservationLawInitialCondition() {
_nbf = 3;
}
dataCacheDouble *newSourceTerm(dataCacheMap &cacheMap)const {
//return new gaussian(cacheMap,0.2,0.3);
return new gaussian(cacheMap,0.5,0.5);
}
};
int main(int argc, char **argv){
GmshMergeFile("input/square1.msh");
//we probably need a class to group those three ones
std::vector<dgGroupOfElements*> elementGroups;
std::vector<dgGroupOfFaces*> faceGroups;
std::vector<dgGroupOfFaces*> boundaryGroups;
int order=1;
int dimension=2;
dgAlgorithm algo;
function::registerDefaultFunctions();
algo.buildGroups(GModel::current(),dimension,order,elementGroups,faceGroups,boundaryGroups);
//for now, we suppose there is only one group of elements
int nbNodes = elementGroups[0]->getNbNodes();
dgConservationLaw *law = dgNewConservationLawWaveEquation();
fullMatrix<double> sol(nbNodes,elementGroups[0]->getNbElements()*law->nbFields());
fullMatrix<double> residu(nbNodes,elementGroups[0]->getNbElements()*law->nbFields());
/*{ // use this block to test functionDerivator
dgConservationLawInitialCondition initLaw;
dataCacheMap cacheMap;
dataCacheElement &cacheElement = cacheMap.getElement();
cacheMap.provideData("UVW").set(elementGroups[0]->getIntegrationPointsMatrix());
cacheElement.set(elementGroups[0]->getElement(0));
dataCacheDouble &xyz = cacheMap.get("XYZ");
dataCacheDouble *source = initLaw.newSourceTerm(cacheMap);
functionDerivator fd (*source, xyz, 1e-12);
xyz().print();
fd.compute();
exit(0);
}*/
// initial condition
{
dgConservationLawInitialCondition initLaw;
algo.residualVolume(initLaw,*elementGroups[0],sol,residu);
algo.multAddInverseMassMatrix(*elementGroups[0],residu,sol);
}
law->addBoundaryCondition("Left",dgNewBoundaryConditionWaveEquationWall());
law->addBoundaryCondition("Right",dgNewBoundaryConditionWaveEquationWall());
law->addBoundaryCondition("Top",dgNewBoundaryConditionWaveEquationWall());
law->addBoundaryCondition("Bottom",dgNewBoundaryConditionWaveEquationWall());
print("output/p.pos",*elementGroups[0],&sol(0,0),0,3);
print_vector("output/uv.pos",*elementGroups[0],&sol(0,0),1,3);
for(int iT=0; iT<10000; iT++) {
algo.rungeKutta(*law,elementGroups,faceGroups,boundaryGroups,1e-3,residu,sol);
if(iT%10==0){
char name[100];
printf("%i\n",iT);
sprintf(name,"output/p_%05i.pos",iT);
print(name,*elementGroups[0],&sol(0,0),0,3);
sprintf(name,"output/uv_%05i.pos",iT);
print_vector(name,*elementGroups[0],&sol(0,0),1,3);
sprintf(name,"output/v_%05i.pos",iT);
print(name,*elementGroups[0],&sol(0,0),2,3);
sprintf(name,"output/u_%05i.pos",iT);
print(name,*elementGroups[0],&sol(0,0),1,3);
}
}
delete law;
}
void print (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields) {
FILE *file = fopen(filename,"w");
fprintf(file,"View \"%s\" {\n", filename);
for(int iel=0;iel<els.getNbElements();iel++){
MElement *el = els.getElement(iel);
fprintf(file,"ST (");
for (int iv=0; iv<el->getNumVertices(); iv++) {
MVertex *vertex = el->getVertex(iv);
SPoint3 coord = vertex->point();
fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
}
fprintf(file,"{");
v+=(iField)*el->getNumVertices();
for (int iv=0; iv<el->getNumVertices(); iv++){
fprintf(file,"%e%c ",*(v++),iv==el->getNumVertices()-1?'}':',');
}
v+=(nbFields-1-iField)*el->getNumVertices();
fprintf(file,";\n");
}
fprintf(file,"};");
fclose(file);
}
void print_vector (const char *filename,const dgGroupOfElements &els, double *v,int iField, int nbFields) {
FILE *file = fopen(filename,"w");
fprintf(file,"View \"%s\" {\n", filename);
for(int iel=0;iel<els.getNbElements();iel++){
MElement *el = els.getElement(iel);
fprintf(file,"VT (");
for (int iv=0; iv<el->getNumVertices(); iv++) {
MVertex *vertex = el->getVertex(iv);
SPoint3 coord = vertex->point();
fprintf(file,"%e, %e, %e%c ",coord.x(),coord.y(),coord.z(),iv==el->getNumVertices()-1?')':',');
}
v+=(iField)*el->getNumVertices();
fprintf(file,"{%e, %e, 0, %e, %e, 0, %e, %e, 0};\n", v[0],v[3],v[1],v[4],v[2],v[5]);
v+=(nbFields-iField)*el->getNumVertices();
}
fprintf(file,"};");
fclose(file);
}
Solver/dgSystemOfEquations.cpp
+ 16
− 3
View file @ b4ffa082
......@@ -60,11 +60,16 @@ dgSystemOfEquations::dgSystemOfEquations(lua_State *L){
// set the conservation law as a string (for now)
int dgSystemOfEquations::setConservationLaw(lua_State *L){
_claw = 0;
//int argc = (int)luaL_checknumber(L,0);
_cLawName = std::string (luaL_checkstring(L, 1));
if (_cLawName == "WaveEquation")
_claw = dgNewConservationLawWaveEquation();
else if (_cLawName == "ShallowWater2d")
_claw = dgNewConservationLawShallowWater2d();
else if (_cLawName == "AdvectionDiffusion"){
std::string advFunction = luaL_checkstring(L,2);
_claw = dgNewConservationLawAdvection(advFunction);
}
if (!_claw)throw;
return 0;
}
......@@ -78,7 +83,15 @@ int dgSystemOfEquations::setOrder(lua_State *L){
int dgSystemOfEquations::addBoundaryCondition (lua_State *L){
std::string physicalName(luaL_checkstring(L, 1));
std::string bcName(luaL_checkstring(L, 2));
if (_cLawName == "WaveEquation"){
//generic boundary conditions
if (bcName == "0Flux"){
_claw->addBoundaryCondition(physicalName,dgBoundaryCondition::new0FluxCondition(*_claw));
}
else if (bcName == "OutsideValues"){
_claw->addBoundaryCondition(physicalName,dgBoundaryCondition::newOutsideValueCondition(*_claw,luaL_checkstring(L,3)));
}
//specific boundary conditions
else if (_cLawName == "WaveEquation"){
if (bcName == "Wall"){
_claw->addBoundaryCondition(physicalName,dgNewBoundaryConditionWaveEquationWall());
}
......@@ -100,6 +113,7 @@ int dgSystemOfEquations::L2Projection (lua_State *L){
_algo->residualVolume(Law,*_elementGroups[i],*_solution->_dataProxys[i],*_rightHandSide->_dataProxys[i]);
_algo->multAddInverseMassMatrix(*_elementGroups[i],*_rightHandSide->_dataProxys[i],*_solution->_dataProxys[i]);
}
return 0;
}
// ok, we can setup the groups and create solution vectors
......@@ -114,7 +128,6 @@ int dgSystemOfEquations::setup(lua_State *L){
// compute the total size of the soution
_solution = new dgDofContainer(_elementGroups,*_claw);
_rightHandSide = new dgDofContainer(_elementGroups,*_claw);
// printf("aaaaaaaaaaaaa\n");
}
......
Solver/luaFunction.cpp
+ 3
− 3
View file @ b4ffa082
......@@ -74,8 +74,7 @@ static int newLuaFunction (lua_State *L){
}
static int newConstantFunction (lua_State *L){
fullMatrix<double> * _ud = (fullMatrix<double> *) lua_touserdata(L, 1);
fullMatrix<double> * _ud = Luna<fullMatrix<double> >::check(L, 1);
int ITER = 1;
std::string functionName = "ConstantFunction";
while (function::get(functionName, true)){
......@@ -91,7 +90,8 @@ static int newConstantFunction (lua_State *L){
int RegisterFunctions (lua_State *L) {
luaL_reg fcts[] = {{"lua", newLuaFunction },
{"constant", newConstantFunction }};
{"constant", newConstantFunction },
{0,0}};
luaL_register(L, "createFunction", fcts);
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment