diff --git a/CMakeLists.txt b/CMakeLists.txt index 8fca86e060a4b410f4d36ac4c51f48e5da47c598..dba191531708f11a42bc2b2dc2d80a862cef508e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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}) diff --git a/Solver/TESTCASES/Advection.lua b/Solver/TESTCASES/Advection.lua new file mode 100644 index 0000000000000000000000000000000000000000..d66da7451f3bf76a2f5302b844f596646b47c26e --- /dev/null +++ b/Solver/TESTCASES/Advection.lua @@ -0,0 +1,45 @@ +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 diff --git a/Solver/TESTCASES/Stommel.lua b/Solver/TESTCASES/Stommel.lua new file mode 100644 index 0000000000000000000000000000000000000000..e420a01be00b27003d001b9b569c2fd2bd6a281d --- /dev/null +++ b/Solver/TESTCASES/Stommel.lua @@ -0,0 +1,18 @@ + +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 diff --git a/Solver/TESTCASES/WavePulse.lua b/Solver/TESTCASES/WavePulse.lua index 74f2c2f70e5c0893e652b5047807d71541ad7955..1e0e93c4a9e3d21bf9fe20cffe062782daef20c8 100644 --- a/Solver/TESTCASES/WavePulse.lua +++ b/Solver/TESTCASES/WavePulse.lua @@ -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 diff --git a/Solver/TESTCASES/stommel_square.geo b/Solver/TESTCASES/stommel_square.geo new file mode 100644 index 0000000000000000000000000000000000000000..56e56036a51d79c2b3b0ab3dbb7bdb599188d22f --- /dev/null +++ b/Solver/TESTCASES/stommel_square.geo @@ -0,0 +1,16 @@ +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}; diff --git a/Solver/dgMainShallowWater2d.cpp b/Solver/dgMainShallowWater2d.cpp deleted file mode 100644 index e3aed3f67da5bcc0bdbff7d9291cf2563ad33023..0000000000000000000000000000000000000000 --- a/Solver/dgMainShallowWater2d.cpp +++ /dev/null @@ -1,71 +0,0 @@ -#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); -} diff --git a/Solver/dgMainTest.cpp b/Solver/dgMainTest.cpp deleted file mode 100644 index bb5c8a25a95ccc877e10076bb5a0fe7af83105be..0000000000000000000000000000000000000000 --- a/Solver/dgMainTest.cpp +++ /dev/null @@ -1,111 +0,0 @@ -#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); -} diff --git a/Solver/dgMainWaveEquation.cpp b/Solver/dgMainWaveEquation.cpp deleted file mode 100644 index 388eba7ee9925a15fbdb617237496a32a78fc8ce..0000000000000000000000000000000000000000 --- a/Solver/dgMainWaveEquation.cpp +++ /dev/null @@ -1,141 +0,0 @@ -#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); -} diff --git a/Solver/dgSystemOfEquations.cpp b/Solver/dgSystemOfEquations.cpp index acaa38827677a29cc3f6f91d85248e3f53f803f3..7f038fa95ebc87241b542ccbaaae16fb450f9e6f 100644 --- a/Solver/dgSystemOfEquations.cpp +++ b/Solver/dgSystemOfEquations.cpp @@ -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") + 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"); } diff --git a/Solver/luaFunction.cpp b/Solver/luaFunction.cpp index 74cde201e6848db93c4dde431028165f71154d1c..0b2f766c655fcdb3c0abb2c49631bddd47055bc4 100644 --- a/Solver/luaFunction.cpp +++ b/Solver/luaFunction.cpp @@ -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); };