diff --git a/ConSol/ConSol/Base/Setup.cpp b/ConSol/ConSol/Base/Setup.cpp index 8e96787..96ebac5 100644 --- a/ConSol/ConSol/Base/Setup.cpp +++ b/ConSol/ConSol/Base/Setup.cpp @@ -1,429 +1,433 @@ // // Setup.cpp // ConSol // // Created by Fabian Moenkeberg on 30.06.17. // Copyright © 2017 EPFL. All rights reserved. // #include "Setup.hpp" Setup::Setup(Configuration config){ dim = config.getDim(); switch (config.getModel()){ case LINEARADVECTION: model = new LinAdv(); break; case BURGERS: model = new Burgers(); break; case BUCKLEYLEVERETT: model = new BuckleyLeverett(); break; case EULER: if (dim == 1){ model = new Euler(); }else{ model = new Euler2D(); } break; case KPP: model = new Kpp(); break; case KPPX: model = new Kppx(); break; case KPPY: model = new Kppy(); break; default: assert(false && "Unknown Model!"); break; } switch (config.getLimiter()){ case NOLIMITER: limiter = new NoLimiter(); break; case EXTREMAPRES: std::cout << "Not implemented limiter" << input_file << std::endl; assert(0 && "Unknown Limiter"); break; default: assert(0 && "Unknown Limiter!"); break; } switch (config.getReconstr()){ case ENO: if (dim ==1){ reconstr = new class ENORec(); }else if (dim == 2){ // reconstr = new class ENO_2D(); switch (config.getStencilChoice()) { case STD: reconstr = new class ENO_2D_std(); break; case OTF: reconstr = new class ENO_2D_otf(); break; case SYMM: reconstr = new class ENO_2D_symm(); break; } } break; case WENO: if (dim ==1){ assert(0 && "Unknown Reconstruction!"); }else if (dim == 2){ if (config.getMeshType() == QUAD) { reconstr = new class WENO_structured(); - }else{ + }else if (config.getMeshType() == HYBRID){ reconstr = new class HybridRBFWENO(); + }else{ + reconstr = new class MRWENO(); } } break; case NORECONSTR: if (dim ==1){ reconstr = new NoReconstr(); }else if (dim == 2) { reconstr = new NoReconstr_2D(); } break; default: assert(0 && "Unknown Reconstruction!"); break; } switch (config.getFlux()){ case LAXFRIEDRICHS: flux = new LaxFried(config.getDim()); break; default: assert(0 && "Unknown Flux!"); break; } // timeInt = new IntBase(config.getTimeIntegration(), flux); switch (config.getTimeIntegration()){ case SSPRK2: timeInt = new Ssprk2(flux); break; case SSPRK3: timeInt = new Ssprk3(flux); break; case FE: timeInt = new Fe(flux); break; case SSPRK5: timeInt = new Ssprk5(flux); break; } switch (config.getInterpBase()){ case MULTIQUADS: if (dim == 1) { interpBase = new Multiquadratics(config.getOrder()-1, config.getShapeParam(), config.getDim()); }else if (dim == 2){ // interpBase = new Multiquadratics_2D(config.getOrder()-1, config.getDegPoly(), config.getShapeParam()); if (config.getMeshType()== MIXED){ interpBase = new Multiquadratics_2D_mixed(config.getOrder()-1, config.getDegPoly(), config.getShapeParam()); }else { switch (config.getStencilChoice()) { case STD: interpBase = new Multiquadratics_2D(config.getOrder() - 1, config.getDegPoly(), config.getShapeParam()); break; case OTF: interpBase = new Multiquadratics_2D_otf(config.getOrder() - 1, config.getDegPoly(), config.getShapeParam()); break; case SYMM: // interpBase = new Multiquadratics_2D_otf(config.getOrder() - 1, config.getDegPoly(), // config.getShapeParam()); - interpBase = new Multiquadratics_2D(config.getOrder() - 1, config.getDegPoly(), - config.getShapeParam()); +// interpBase = new Multiquadratics_2D(config.getOrder() - 1, config.getDegPoly(), +// config.getShapeParam()); + interpBase = new Multiquadratics_2D_WENO(config.getOrder() - 1, config.getDegPoly(), + config.getShapeParam()); break; } } } break; case POLYHARMONICS: if (dim == 1) { interpBase = new Polyharmonics(config.getOrder()-1, config.getDim()); }else if (dim == 2){ interpBase = new Polyharmonics_2D(config.getOrder()-1); } break; case GAUSSIANS: if (dim == 1) { assert(0 && "Gaussians exist only for 2 space dimensions!"); }else if (dim == 2){ interpBase = new Gaussians_2D(config.getOrder()-1); } break; default: assert(0 && "Unknown Interpolation Base!"); break; } // Construct Mesh switch (config.getBc()){ case PERIODIC: if (dim==1){ bc = new Periodic();} else bc = new Periodic2(); break; case ABSORBINGBC: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new AbsorbingBC(); break; case MIXEDBC: if (dim==1){ assert(false && "Not Implemented, yet!"); bc = new MixedBC();} else bc = new MixedBC_2D(); break; case NEUMANN: assert(false && "Not Implemented, yet!"); break; case SLIPWALL: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new SlipWallBC(); break; case DOUBLEMACHREF: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new DoubleMachReflectionBC(); break; case MIXEDBCCOMPONENT: assert(false && "Not Implemented, yet!"); break; case DIRICHLET: assert(false && "Not Implemented, yet!"); break; case FAC: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new BCFAC(); break; case OPENFLOW: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new OpenFlowBC(); break; case SCRAMJET: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new ScramJetBC(); break; } switch (config.getDim()){ case 1:{ double delta = (config.getXmax()-config.getXmin())/config.getNx()[0]; std::vector Nx = config.getNx(); std::vector x0(Nx[0]+1); x0.resize(Nx[0]+2); for (int i = 0; i < Nx[0]+2; i++){ x0[i] = i*delta - 0.5/Nx[0]; } int nX0 = (int)x0.size(); std::vector> x(nX0-1,std::vector(1)); for (int i = 0; i < nX0-1; i++){ x[i][0] = x0[i+1]; } mesh = new Rect(x); break;} case 2:{ if (!config.getLoadGrid()){ assert(false && "Two Dim grids must get loaded from file."); } switch (config.getMeshType()) { case TRI: mesh = new Tri(); break; case HYBRID: mesh = new FastHybridMesh(); break; case QUAD: mesh = new StructuredMesh(); break; case MIXED: mesh = new HybridMesh(); break; default: assert(0 && "Unknown Mesh Type!"); break; } break;} } switch (config.getSource()){ case NOSOURCE: source = new Source(); break; } switch (config.getInitialCond()){ case SIN: initialCond = new Sin(); break; case BUBBLE: initialCond = new bubble(); break; case COSCOS: initialCond = new coscos(); break; case SINSIN: initialCond = new sinsin(); break; case InitKPP: initialCond = new initKPP(); break; case InitKPPY: initialCond = new initKPPY(); break; case EULERSOD: initialCond = new InitEulerSod(); break; case InitBURGERS: initialCond = new initBurgers(); break; case SHOCKX: initialCond = new shockX(); break; case InitCUBE: initialCond = new cube(); break; case InitISENTROPICVORTEX: initialCond = new initIsentropicVortex(); break; case InitEULERBLAST: initialCond = new initEulerBlast(); break; case InitVORTEXSHOCK: initialCond = new initVortexShock(); break; case RP: initialCond = new initRP(config.getConstants()); break; case InitDOUBLEMACHREF: initialCond = new initDoubleMachReflection(); break; case InitFAC: if (config.getBc() == SCRAMJET){ initialCond = new initFAC(0, 3.0); }else { initialCond = new initFAC(); } break; case InitAIRFOIL: initialCond = new AirFoil(); break; } } Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Source *source, Reconstruct *reconstr, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){ this->mesh = mesh; this->model = model; this->timeInt = timeInt; this->source = source; this->reconstr = reconstr; this->flux = flux; this->initialCond = initialCond; this->limiter = limiter; this->interpBase = new InterpBase(); this->bc = bc; } Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Source *source, Reconstruct *reconstr, InterpBase *interpBase, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){ this->mesh = mesh; this->model = model; this->timeInt = timeInt; this->source = source; this->reconstr = reconstr; this->flux = flux; this->initialCond = initialCond; this->limiter = limiter; this->interpBase = interpBase; this->bc = bc; } Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Reconstruct *reconstr, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){ this->mesh = mesh; this->model = model; this->timeInt = timeInt; this->source = new Source(); this->reconstr = reconstr; this->flux = flux; this->initialCond = initialCond; this->limiter = limiter; this->interpBase = new InterpBase(); this->bc = bc; } Setup::Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Reconstruct *reconstr, InterpBase *interpBase, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond){ this->mesh = mesh; this->model = model; this->timeInt = timeInt; this->source = new Source(); this->reconstr = reconstr; this->flux = flux; this->initialCond = initialCond; this->limiter = limiter; this->interpBase = interpBase; this->bc = bc; } Solution* Setup::initialize(Configuration config){ model->initialize(config); mesh->initialize(bc, config); interpBase->initialize(config); reconstr->initialize(model, interpBase, limiter, mesh, config); flux->initialize(mesh, source, model, reconstr); limiter->initialize(interpBase, model); if (source->getIsUsed()){ source->initialize(config); } timeInt->initialize(config, mesh); this->soln = Solution(config, mesh, initialCond); return &soln; } MeshBase* Setup::getMesh(){return mesh;}; BcBase* Setup::getBC() {return bc;}; NumFlux* Setup::getFlux(){return flux;}; ModelBase* Setup::getModel(){return model;}; IntBase* Setup::getTimeInt(){return timeInt;}; LimiterBase* Setup::getLimiter(){return limiter;}; void Setup::exportStencils(std::vector> &U, MeshBase *mesh, double t, Configuration config){ reconstr->exportStencils(U, mesh, t, config.getOutputPath()); } diff --git a/ConSol/ConSol/Base/Setup.hpp b/ConSol/ConSol/Base/Setup.hpp index 0dcfc75..99427bf 100644 --- a/ConSol/ConSol/Base/Setup.hpp +++ b/ConSol/ConSol/Base/Setup.hpp @@ -1,202 +1,204 @@ // // Setup.hpp // ConSol // // Created by Fabian Moenkeberg on 30.06.17. // Copyright © 2017 EPFL. All rights reserved. // #ifndef Setup_hpp #define Setup_hpp #include #include "../Mesh/MeshBase.hpp" #include "../Model/ModelBase.hpp" #include "../TimeIntegration/IntBase.hpp" #include "../TimeIntegration/Ssprk2.hpp" #include "../TimeIntegration/Ssprk3.hpp" #include "../TimeIntegration/Ssprk5.hpp" #include "../TimeIntegration/Fe.hpp" #include "../Source/Source.hpp" #include "../Reconstr/Reconstruct.hpp" #include "../Flux/NumFlux.hpp" #include "Solution.hpp" #include "../Utility/AbstractFunct.hpp" #include "../Reconstr/Limiter/LimiterBase.hpp" #include "../Reconstr/InterpolationBase/InterpBase.hpp" #include "../Mesh/BC/BcBase.hpp" #include #include "../Reconstr/Limiter/NoLimiter.hpp" #include "../Reconstr/Limiter/ExtremaPresLimiter.hpp" #include "../Reconstr/NoReconstr.hpp" #include "../Reconstr/NoReconstr_2D.hpp" #include "../Mesh/Rect.hpp" #include "../Mesh/Tri.h" #include "../Mesh/HybridMesh.hpp" #include "../Mesh/FastHybridMesh.hpp" #include "../Mesh/StructuredMesh.hpp" #include "../Flux/LaxFried.hpp" #include "../Flux/LaxFriedGlobal.hpp" #include "../Model/LinAdv.hpp" #include "../Model/Burgers.hpp" #include "../Model/Euler.hpp" #include "../Model/Euler2D.hpp" #include "../Model/BuckleyLeverett.hpp" #include "../Model/Kpp.hpp" #include "../Model/Kppx.hpp" #include "../Model/Kppy.hpp" #include "../Utility/InitialConditions/Sin.hpp" #include "../Utility/InitialConditions/2D/bubble.hpp" #include "../Utility/InitialConditions/2D/InitKPP.hpp" #include "../Utility/InitialConditions/2D/InitKPPY.hpp" #include "../Utility/InitialConditions/InitEulerSod.hpp" #include "../Utility/InitialConditions/2D/coscos.hpp" #include "../Utility/InitialConditions/2D/sinsin.hpp" #include "../Utility/InitialConditions/2D/InitBurgers.hpp" #include "../Utility/InitialConditions/2D/shockX.hpp" #include "../Utility/InitialConditions/2D/cube.hpp" #include "../Utility/InitialConditions/2D/InitIsentropicVortex.hpp" #include "../Utility/InitialConditions/2D/InitEulerBlast.hpp" #include "../Utility/InitialConditions/2D/InitVortexShock.hpp" #include "../Utility/InitialConditions/2D/InitFAC.hpp" #include "../Utility/InitialConditions/2D/AirFoil.hpp" #include "../Utility/InitialConditions/2D/RP.hpp" #include "../Utility/InitialConditions/2D/InitDoubleMachReflection.hpp" #include "../Reconstr/InterpolationBase/Multiquadratics.hpp" #include "../Reconstr/InterpolationBase/Multiquadratics_2D.hpp" #include "../Reconstr/InterpolationBase/Multiquadratics_2D_mixed.hpp" #include "../Reconstr/InterpolationBase/Multiquadratics_2D_otf.hpp" +#include "../Reconstr/InterpolationBase/Multiquadratics_2D_WENO.hpp" #include "../Reconstr/InterpolationBase/Polyharmonics.hpp" #include "../Reconstr/InterpolationBase/Polyharmonics_2D.hpp" #include "../Reconstr/InterpolationBase/Gaussians_2D.hpp" #include "../Reconstr/ENO.hpp" #include "../Reconstr/ENO_2D.hpp" #include "../Reconstr/ENO_2D_otf.hpp" #include "../Reconstr/ENO_2D_std.hpp" #include "../Reconstr/ENO_2D_symm.hpp" +#include "../Reconstr/MRWENO.hpp" #include "../Reconstr/CWENO_2D.h" #include "../Reconstr/WENO_structured.hpp" #include "../Reconstr/HybridRBFWENO.hpp" #include "../Mesh/BC/Periodic.hpp" #include "../Mesh/BC/Periodic2.hpp" #include "../Mesh/BC/AbsorbingBC.hpp" #include "../Mesh/BC/MixedBC.hpp" #include "../Mesh/BC/MixedBC_2D.hpp" #include "../Mesh/BC/SlipWallBC.hpp" #include "../Mesh/BC/BCFAC.hpp" #include "../Mesh/BC/DoubleMachReflectionBC.hpp" #include "../Mesh/BC/OpenFlowBC.hpp" #include "../Mesh/BC/ScramJetBC.hpp" class Setup{ MeshBase *mesh; ModelBase *model; IntBase *timeInt; Source *source; Reconstruct *reconstr; NumFlux *flux; Solution soln; AbstractFunct *initialCond; LimiterBase *limiter; InterpBase *interpBase; BcBase *bc; ModelName modelName; evaluationMethod evalMethod; bcName bcname; bcName bcType[6]; ReconstrName reconstrName; InterpBaseName interpBaseName; TimeIntegrationName time_scheme; FluxName fluxScheme; limiterName limiter_name; bool reconstrOnCharact = 0; int order; double tMax; double cfl; double shapeParam; int maxNumWrite = 100; int max_iter = 50000; double min_residue; std::vector tInclude; double xmin,xmax,ymin,ymax; int dim; std::vector Nx_cell, Ny_cell; char* input_file; std::string ic_file; std::string source_file; std::string OUTPUT_PATH; std::string SOLN_DIR; bool ic_file_read; bool exact_available; public: Setup(){ }; ~Setup(){ delete mesh; delete model; delete timeInt; delete flux; delete reconstr; delete limiter; delete source; delete interpBase; delete bc; delete initialCond; // delete interpBase; }; // Setup(char* file); Setup(Configuration config); Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Source *source, Reconstruct *reconstr, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond); Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Source *source, Reconstruct *reconstr, InterpBase *interpBase, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond); Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Reconstruct *reconstr, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond); Setup(MeshBase *mesh, BcBase *bc, ModelBase *model, IntBase *timeInt, Reconstruct *reconstr, InterpBase *interpBase, LimiterBase *limiter, NumFlux *flux, AbstractFunct *initialCond); Solution* initialize(Configuration config); MeshBase* getMesh(); BcBase* getBC(); NumFlux* getFlux(); ModelBase* getModel(); IntBase* getTimeInt(); LimiterBase* getLimiter(); void exportStencils(std::vector> &U, MeshBase *mesh, double t, Configuration config); // void read (); // Configuration* getConfig(); private: // void read_grid (Reader&); // void read_constants (Reader&); // void read_numeric (Reader&); // void read_model (Reader&); // void read_initial_condition (Reader&); // void read_exact (Reader&); // void read_boundary_condition (Reader&); // void read_output (Reader&); void process_param(Configuration config); }; #endif /* Setup_hpp */ diff --git a/ConSol/ConSol/CMakeLists.txt b/ConSol/ConSol/CMakeLists.txt index 2a5e9ff..5f1ae87 100644 --- a/ConSol/ConSol/CMakeLists.txt +++ b/ConSol/ConSol/CMakeLists.txt @@ -1,286 +1,286 @@ find_package(Threads REQUIRED) if (APPLE) cmake_policy(SET CMP0074 NEW) set(HDF5_ROOT "/usr/local/Cellar/hdf5/1.10.5") ## set(ENV{GMSH_INC} "/usr/local/Cellar/gmsh/4.4.1/share/doc/gmsh/demos/api") ## set(ENV{GMSH_LIB} "/usr/local/Cellar/gmsh/4.4.1/lib") else() ### NOTE: on moenkerg-laptop, the next two variables should NOT be set!!! set(ENV{HDF5_ROOT} "/home/moenkebe/lib_loc/hdf5") set(Eigen3_DIR "/home/moenkebe/lib_loc/eigen3") set(ENV{OpenBLAS_DIR} "/home/moenkebe/lib_loc/openblas") # set(ENV{GMSH_INC} "/home/moenkebe/lib_loc/gmsh/include") # set(ENV{GMSH_LIB} "/home/moenkebe/lib_loc/gmsh/lib") # set(ENV{GMSH_DIR} "/home/moenkebe/lib_loc/gmsh") endif() set(HDF5_VALID_LANGUAGE_BINDINGS CXX) find_package(HDF5 REQUIRED COMPONENTS CXX) include_directories(${HDF5_INCLUDE_DIRS}) ### Find GMSH if (APPLE) find_library(GMSH_LIB gmsh) if(NOT GMSH_LIB) message(FATAL_ERROR "Could not find libgmsh") endif() find_path(GMSH_INC gmsh.h) if(NOT GMSH_INC) message(FATAL_ERROR "Could not find gmsh.h") endif() else() find_library(GMSH_LIB gmsh HINTS "/home/moenkebe/lib_loc/gmsh/lib") if(NOT GMSH_LIB) message(FATAL_ERROR "Could not find libgmsh") endif() find_path(GMSH_INC gmsh.h "/home/moenkebe/lib_loc/gmsh/include") if(NOT GMSH_INC) message(FATAL_ERROR "Could not find gmsh.h") endif() endif() if(GMSH_LIB MATCHES ".a") # FIXME - generalize this find_library(BLAS_LIB blas) if(BLAS_LIB) list(APPEND EXTRA_LIBS ${BLAS_LIB}) endif() find_library(LAPACK_LIB lapack) if(LAPACK_LIB) list(APPEND EXTRA_LIBS ${LAPACK_LIB}) endif() endif() include_directories(${GMSH_INC}) ### NOTE: on moenkerg-laptop, the next two lines should be commented!!! set(USE_OPENMP 1) set(OPENBLAS_NUM_THREADS 1) set(GOTO_NUM_THREADS 1) set(OMP_NUM_THREADS 1) message(${USE_OPENMP}) find_package( OpenBLAS REQUIRED) include_directories(${OpenBLAS_INCLUDE_DIRS}) if (UNIX AND NOT APPLE) ### NOTE: next two lines need to be commented on mathicsepc13 and 55 # find_package( LAPACKE REQUIRED ) # include_directories(${LAPACKE_INCLUDE_DIR}) endif() find_package (Eigen3 3.3 REQUIRED NO_MODULE) add_library(UtilityLib STATIC Utility/triangle_dunavant_rule.cpp Utility/triangle_dunavant_rule.hpp Utility/gauss-legendre-quadrature/gauss_legendre.cpp Utility/gauss-legendre-quadrature/gauss_legendre.hpp Utility/JBurkardt/triangle_integrals.cpp Utility/JBurkardt/triangle_integrals.hpp Utility/Quadrature1D/quadrule.cpp Utility/Quadrature1D/quadrule.hpp Utility/JBurkardt/triangle_wandzura_rule.cpp Utility/JBurkardt/triangle_wandzura_rule.hpp Utility/JBurkardt/condition.cpp Utility/JBurkardt/condition.hpp Utility/JBurkardt/condition_prb.cpp Utility/JBurkardt/r8lib.cpp Utility/JBurkardt/r8lib.hpp) add_library(Model Base/Configuration.cpp Base/Configuration.hpp Base/runSolver.cpp Base/runSolver.hpp Base/Setup.cpp Base/Setup.hpp Base/SolnMap.cpp Base/SolnMap.hpp Base/Solution.cpp Base/Solution.hpp Base/reader.h Flux/LaxFried.cpp Flux/LaxFried.hpp Flux/LaxFriedGlobal.cpp Flux/LaxFriedGlobal.hpp Flux/NetFlux.cpp Flux/NetFlux.hpp Flux/NumFlux.cpp Flux/NumFlux.hpp Mesh/BC/BcBase.cpp Mesh/BC/BcBase.hpp Mesh/BC/MixedBC.cpp Mesh/BC/MixedBC.hpp Mesh/BC/Periodic.cpp Mesh/BC/Periodic.hpp Mesh/MeshBase.cpp Mesh/MeshBase.hpp Mesh/Rect.cpp Mesh/Rect.hpp Model/BuckleyLeverett.cpp Model/BuckleyLeverett.hpp Model/Burgers.cpp Model/Burgers.hpp Model/Euler.cpp Model/Euler.hpp Model/LinAdv.cpp Model/LinAdv.hpp Model/ModelBase.cpp Model/ModelBase.hpp Reconstr/InterpolationBase/InterpBase.cpp Reconstr/InterpolationBase/InterpBase.hpp Reconstr/InterpolationBase/Multiquadratics.cpp Reconstr/InterpolationBase/Multiquadratics.hpp Reconstr/InterpolationBase/Polyharmonics.cpp Reconstr/InterpolationBase/Polyharmonics.hpp Reconstr/Limiter/ExtremaPresLimiter.cpp Reconstr/Limiter/ExtremaPresLimiter.hpp Reconstr/Limiter/LimiterBase.cpp Reconstr/Limiter/LimiterBase.hpp Reconstr/Limiter/NoLimiter.cpp Reconstr/Limiter/NoLimiter.hpp Reconstr/ENO.cpp Reconstr/ENO.hpp Reconstr/NoReconstr.cpp Reconstr/NoReconstr.hpp Reconstr/Reconstruct.cpp Reconstr/Reconstruct.hpp Source/Source.cpp Source/Source.hpp TimeIntegration/IntBase.cpp TimeIntegration/IntBase.hpp Utility/InitialConditions/InitBuckleyLeverett.cpp Utility/InitialConditions/InitBuckleyLeverett.hpp Utility/InitialConditions/InitEulerLaxShockTube.cpp Utility/InitialConditions/InitEulerLaxShockTube.hpp Utility/InitialConditions/InitEulerLaxShockTube2.cpp Utility/InitialConditions/InitEulerLaxShockTube2.hpp Utility/InitialConditions/InitEulerShockEntropy.cpp Utility/InitialConditions/InitEulerShockEntropy.hpp Utility/InitialConditions/InitEulerSod.cpp Utility/InitialConditions/InitEulerSod.hpp Utility/InitialConditions/InitEulerWCBlast.cpp Utility/InitialConditions/InitEulerWCBlast.hpp Utility/InitialConditions/Sin.cpp Utility/InitialConditions/Sin.hpp Utility/AbstractFunct.cpp Utility/AbstractFunct.hpp Utility/calcError.cpp Utility/calcError.hpp Utility/ddGeneralizedNewton.cpp Utility/ddGeneralizedNewton.hpp Utility/findNearest.cpp Utility/findNearest.hpp Utility/GetQuadraturePoints.cpp Utility/GetQuadraturePoints.hpp Utility/isWellDef.cpp Utility/isWellDef.hpp Utility/MatrixVectorProd.cpp Utility/MatrixVectorProd.hpp Utility/signPropFulfilled.cpp Utility/signPropFulfilled.hpp Utility/vvra.cpp Utility/vvra.hpp Utility/vvraEigen.cpp Utility/vvraEigen.hpp Base/reader.h Utility/process_command_line.cpp Utility/progress_display.cpp Mesh/Tri.cpp Mesh/Tri.h Mesh/gmsh_io.cpp Mesh/gmsh_io.hpp Mesh/BC/Periodic2.cpp Mesh/BC/Periodic2.hpp Flux/LaxFlux.cpp Flux/LaxFlux.hpp Flux/LaxFlux_2D.cpp Flux/LaxFlux_2D.hpp Utility/InitialConditions/2D/bubble.cpp Utility/InitialConditions/2D/bubble.hpp Reconstr/NoReconstr_2D.cpp Reconstr/NoReconstr_2D.hpp Reconstr/ENO_2D.cpp Reconstr/ENO_2D.hpp Reconstr/InterpolationBase/Multiquadratics_2D.cpp Reconstr/InterpolationBase/Multiquadratics_2D.hpp Model/Kpp.cpp Model/Kpp.hpp Utility/InitialConditions/2D/InitKPP.cpp Utility/InitialConditions/2D/InitKPP.hpp Mesh/BC/MixedBC_2D.cpp Mesh/BC/MixedBC_2D.hpp TimeIntegration/Ssprk3.cpp TimeIntegration/Ssprk3.hpp TimeIntegration/Fe.cpp TimeIntegration/Fe.hpp TimeIntegration/Ssprk2.cpp TimeIntegration/Ssprk2.hpp Utility/InitialConditions/2D/coscos.cpp Utility/InitialConditions/2D/coscos.hpp Model/Kppx.cpp Model/Kppx.hpp Model/Kppy.cpp Model/Kppy.hpp Utility/InitialConditions/2D/InitKPPY.cpp Utility/InitialConditions/2D/InitKPPY.hpp Reconstr/CWENO_2D.cpp Reconstr/CWENO_2D.h Reconstr/ENO_2D_otf.cpp Reconstr/ENO_2D_otf.hpp Reconstr/ENO_2D_std.cpp Reconstr/ENO_2D_std.hpp Utility/InitialConditions/2D/InitBurgers.cpp Utility/InitialConditions/2D/InitBurgers.hpp Utility/InitialConditions/2D/shockX.cpp Utility/InitialConditions/2D/shockX.hpp TimeIntegration/Ssprk5.cpp TimeIntegration/Ssprk5.hpp Utility/InitialConditions/2D/sinsin.cpp Utility/InitialConditions/2D/sinsin.hpp Model/Euler2D.cpp Model/Euler2D.hpp Utility/InitialConditions/2D/InitIsentropicVortex.cpp Utility/InitialConditions/2D/InitIsentropicVortex.hpp Utility/InitialConditions/2D/InitEulerBlast.cpp Utility/InitialConditions/2D/InitEulerBlast.hpp Utility/InitialConditions/2D/InitVortexShock.cpp Utility/InitialConditions/2D/InitVortexShock.hpp Utility/InitialConditions/2D/RP.cpp Utility/InitialConditions/2D/RP.hpp Reconstr/InterpolationBase/Polyharmonics_2D.cpp Reconstr/InterpolationBase/Polyharmonics_2D.hpp Reconstr/InterpolationBase/Multiquadratics_2D_otf.cpp Reconstr/InterpolationBase/Multiquadratics_2D_otf.hpp Mesh/BC/SlipWallBC.cpp Mesh/BC/SlipWallBC.hpp Utility/InitialConditions/2D/InitDoubleMachReflection.cpp Utility/InitialConditions/2D/InitDoubleMachReflection.hpp Mesh/BC/DoubleMachReflectionBC.cpp Mesh/BC/DoubleMachReflectionBC.hpp Mesh/BC/flux_farfield.cpp Mesh/BC/flux_farfield.hpp Reconstr/ENO_2D_symm.cpp Reconstr/ENO_2D_symm.hpp Reconstr/ENO_2D_otf_flex.cpp Reconstr/ENO_2D_otf_flex.hpp Reconstr/InterpolationBase/Gaussians_2D.cpp Reconstr/InterpolationBase/Gaussians_2D.hpp Mesh/HybridMesh.cpp Mesh/HybridMesh.hpp Mesh/StructuredMesh.cpp Mesh/StructuredMesh.hpp Reconstr/WENO_structured.cpp Reconstr/WENO_structured.hpp Utility/Hesthaven/stdWenoWeights.cpp - Utility/Hesthaven/stdWenoWeights.hpp Mesh/FastHybridMesh.cpp Mesh/FastHybridMesh.hpp Reconstr/HybridRBFWENO.cpp Reconstr/HybridRBFWENO.hpp Mesh/BC/AbsorbingBC.cpp Mesh/BC/AbsorbingBC.hpp Utility/InitialConditions/2D/cube.cpp Utility/InitialConditions/2D/cube.hpp Reconstr/InterpolationBase/Multiquadratics_2D_mixed.cpp Reconstr/InterpolationBase/Multiquadratics_2D_mixed.hpp Utility/InitialConditions/2D/InitFAC.cpp Utility/InitialConditions/2D/InitFAC.hpp Mesh/BC/BCFAC.cpp Mesh/BC/BCFAC.hpp Utility/InitialConditions/2D/AirFoil.cpp Utility/InitialConditions/2D/AirFoil.hpp Mesh/BC/OpenFlowBC.cpp Mesh/BC/OpenFlowBC.hpp Mesh/BC/ScramJetBC.cpp Mesh/BC/ScramJetBC.hpp) + Utility/Hesthaven/stdWenoWeights.hpp Mesh/FastHybridMesh.cpp Mesh/FastHybridMesh.hpp Reconstr/HybridRBFWENO.cpp Reconstr/HybridRBFWENO.hpp Mesh/BC/AbsorbingBC.cpp Mesh/BC/AbsorbingBC.hpp Utility/InitialConditions/2D/cube.cpp Utility/InitialConditions/2D/cube.hpp Reconstr/InterpolationBase/Multiquadratics_2D_mixed.cpp Reconstr/InterpolationBase/Multiquadratics_2D_mixed.hpp Utility/InitialConditions/2D/InitFAC.cpp Utility/InitialConditions/2D/InitFAC.hpp Mesh/BC/BCFAC.cpp Mesh/BC/BCFAC.hpp Utility/InitialConditions/2D/AirFoil.cpp Utility/InitialConditions/2D/AirFoil.hpp Mesh/BC/OpenFlowBC.cpp Mesh/BC/OpenFlowBC.hpp Mesh/BC/ScramJetBC.cpp Mesh/BC/ScramJetBC.hpp Reconstr/MRWENO.cpp Reconstr/MRWENO.hpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO.hpp) if (APPLE) target_link_libraries(Model ${OpenBLAS_LIBRARIES} Eigen3::Eigen ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES} UtilityLib ${GMSH_LIB}) else() target_link_libraries(Model ${OpenBLAS_LIBRARIES} ${LAPACKE_LIBRARIES} ${BLAS_LIBRARIES} ${LAPACK_LIBRARIES} Eigen3::Eigen ${HDF5_CXX_LIBRARIES} ${HDF5_LIBRARIES} UtilityLib ${GMSH_LIB}) endif() diff --git a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp new file mode 100644 index 0000000..b98c02c --- /dev/null +++ b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp @@ -0,0 +1,256 @@ +// +// Created by Fabian Moenkeberg on 2019-12-16. +// + +#include "Multiquadratics_2D_WENO.hpp" +#include "cblas.h" + +Multiquadratics_2D_WENO::Multiquadratics_2D_WENO(int order, double eps){ + this->order = order; + this->order_scheme = 2*((order+1)/2 + (order+1)%2); + this->shapeParam = eps; + this->dim = 2; + this->degPoly = order -1; + this->name = MULTIQUADS; + setMulti_Index(degPoly+1); +} + +Multiquadratics_2D_WENO::Multiquadratics_2D_WENO(int order, int degPoly, double eps){ + this->order = 1; +// this->order_scheme = 2*((order+1)/2 + (order+1)%2); + this->order_scheme = order+1; + this->shapeParam = eps; + this->dim = 2; + this->degPoly = degPoly; + this->name = MULTIQUADS; + setMulti_Index(degPoly+1); +} + +std::vector Multiquadratics_2D_WENO::evaluate(std::vector> evPoints, MeshBase *mesh, std::vector stencil, std::vector values){ + int nStencils = (int) stencil.size(); + std::vector>> xB(nStencils,std::vector>(NODE_NUM, std::vector(2))); + std::vector> xi(NODE_NUM,std::vector(2)); + for (int i = 0; i < nStencils; i++){ + xi = mesh->getXCoord(mesh->getCells(stencil[i])); + xB[i] = xi; + } + + double rad = stableRad / mesh->getApproxStencilSize(stencil); +// double rad = stableRad / sqrt(mesh->getSize(cellIndex)); + + std::vector ret = Multiquadratics_2D_WENO::evaluateDirect(evPoints, xB, values, 1, rad); + return ret; +} + +std::vector Multiquadratics_2D_WENO::evaluateDirect(std::vector> evPoints, + std::vector>> xB, + std::vector values, double eps, double rad){ + double eps0 = 1e-10; + int nEv = (int) evPoints.size(); + std::vector smInd(3); + std::vector> ret(3, std::vector(nEv,0)); + int nRbf = (int) xB.size(); + std::vector coeff_old; + int nRbf_old; + int nP_old; + std::vector>> xB_; +// std::vector>> xB_.reserve(nRbf,std::vector>(3, std::vector(2))); + std::vector values_; + std::vector values_old; + int degPolyMax_ = std::max(getDegPoly(nRbf),0); + for (int degPoly_ = 0; degPoly_ <= degPolyMax_; degPoly_++) { + int nP = multi_index.size(); + int nP_tmp = 0; + int sum_multi_index = multi_index[0][0] + multi_index[0][1]; + for (int i = 0; i < nP; i++) { + if (multi_index[i][0] + multi_index[i][1] <= degPoly_) { + nP_tmp++; + } + } + nP = nP_tmp; + + if (degPoly_ == 0) { + nRbf = 1; + xB_.resize(nRbf,std::vector>(3, std::vector(2))); + xB_[0] = xB[0]; + values_.resize(nRbf); + values_[0] = values[0]; + }else if (degPoly_ == 1) { + nRbf = 3; + xB_.resize(nRbf,std::vector>(3, std::vector(2))); + values_.resize(nRbf); + for(int i = 0; i < nRbf; i++){ + xB_[i] = xB[i]; + values_[i] = values[i]; + } + }else{ + nRbf = xB.size(); + xB_.resize(nRbf,std::vector>(3, std::vector(2))); + values_.resize(nRbf); + for(int i = 0; i < nRbf; i++){ + xB_[i] = xB[i]; + values_[i] = values[i]; + } + } + + int m = nRbf + nP; + eps = eps * rad; + + int approxOrder = std::min(degPoly_ + 1, order); + getInterpCoeff(xB_, values_, eps, nP, approxOrder); + std::vector values_tmp = values_; + if (nRbf >1){ + for ( int i = 0; i < nRbf_old; i++){ + values_tmp[i] -= values_old[i]; + } + + for ( int i = 0; i < nP_old; i++){ + values_tmp[i + nRbf] -= values_old[i + nRbf_old]; + } + } + values_old = values_; + std::vector> EvA0(nEv, std::vector(nRbf)); + MVEvaluationMatrix(evPoints, xB_, eps, EvA0, approxOrder); + + std::vector> PA(nEv, std::vector(nP)); + + MVPolynomialEvaluationMatrixV(evPoints, xB_, PA, eps, nP); + + double *EvA = new double[nEv * m]; + + for (int i = 0; i < nEv; i++) { + for (int j = 0; j < nRbf; j++) { + EvA[i + j * nEv] = EvA0[i][j]; + } + for (int j = 0; j < nP; j++) { + EvA[i + (j + nRbf) * nEv] = PA[i][j]; + } + } + + if (nRbf > 1){ + smInd[degPoly_] = 1.0/(ptwiseSmInd(xB[0], xB_, values_tmp, eps) + eps0) + 1; + }else{ +// double a1,a2,a3,b1,b2,b3; +// if (degPolyMax_>0){ +// for (int i = 0; i < 3; i++){ +// +// } +// } + smInd[0] = 1.0/eps0+1; + } + double *sVecPt = values_.data(); + double *y = new double[nEv]; + cblas_dgemv(CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, (int) nEv, (int) (nRbf + nP), 1.0, EvA, (int) nEv, + sVecPt, 1, 0.0, y, (int) 1); + + + for (int i = 0; i < nEv; i++){ + ret[degPoly_][i] = y[i]; + } + +// if (nRbf == 6){ +// getInterpCoeff(xB_, values, eps, nP, approxOrder); +// double *sVecPt = values.data(); +// cblas_dgemv(CblasColMajor, CBLAS_TRANSPOSE::CblasNoTrans, (int) nEv, (int) (nRbf + nP), 1.0, EvA, (int) nEv, +// sVecPt, 1, 0.0, y, (int) 1); +// +// for (int i = 0; i < nEv; i++){ +// for ( int j = 1; j < degPolyMax_; j++) { +// ret[0][i] += ret[j][i]; +// } +// ret[1][i] = y[i]; +// } +// } + +// values_old = values_; + nP_old = nP; + nRbf_old = nRbf; + + delete[] y; + + delete[] EvA; +// delete[] sVecPt; + } + double sum_smInd = 0; + for ( int i = 0; i <= degPolyMax_; i++){ + sum_smInd += smInd[i]; + } + + for ( int j = 0; j < degPolyMax_; j++){ + for (int i = 0; i < nEv; i++){ + ret[degPolyMax_-j][i] += (ret[degPolyMax_-j][i] - ret[degPolyMax_-j-1][i])*smInd[degPolyMax_-j]/sum_smInd; +// ret[degPolyMax_-j][i] += ret[degPolyMax_-j][i]*smInd[degPolyMax_-j]/sum_smInd; + } + } + + for (int i = 0; i < nEv; i++){ + ret[0][i] = ret[0][i]*smInd[0]/sum_smInd; + } + + std::cout << smInd[0]/sum_smInd << ", " << smInd[1]/sum_smInd << ", " << smInd[2]/sum_smInd << "\n"; + return ret[0]; +} + +double Multiquadratics_2D_WENO::ptwiseSmInd(std::vector> xBi, + std::vector>> xB, std::vector values, + double eps){ + + int nRbf = (int) xB.size(); + int degPoly_ = std::max(getDegPoly(nRbf),0); + + double Ind = 0; +// double s = 0; +// double p = 1; +// double dist = 0; + std::vector xM(2, 0.0); + + xM[0] = (xBi[0][0] + xBi[1][0] + xBi[2][0])/NODE_NUM; + xM[1] = (xBi[0][1] + xBi[1][1] + xBi[2][1])/NODE_NUM; + +// double delta; +// double d_0, d_1; + for (int i = 0; i < nRbf;i++){ +// Ind += std::pow(values[i],2); +// delta = eps*std::sqrt(calcSize(xB[i])); + Ind += values[i]*values[i]; +// p *= delta*delta; +// s += delta*delta; +// d_0 = xM[0]-calcMidPt(xB[i],0); +// d_1 = xM[1]-calcMidPt(xB[i],1); +// dist += eps*std::sqrt(d_0*d_0 + d_1*d_1); + } +// Ind = Ind*p*p/(s*s); +// Ind = 0; + Ind/=nRbf; + double delta0 = std::sqrt(calcSize(xBi)); + for (int j = 1;j <= degPoly_; j++){ +// double Ds = 0; +// for (int i = 0; i < nRbf;i++){ +// Ds += values[i]*DF(std::vector{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j); +// } +// Ds += DPolynomial(xM, xB, values, j, values.size() - nRbf, eps); +// Ind += std::pow(Ds,2); + for (int j2 = 0; j2 <= j; j2++){ + double Ds = 0; +// for (int i = 0; i < nRbf;i++){ +// Ds += values[i]*DF(std::vector{eps*(xM[0]-calcMidPt(xB[i],0)),eps*(xM[1]-calcMidPt(xB[i],1))},j2, j-j2, eps); +// } + Ds += DPolynomial(xM, xB, values, j2, j-j2, values.size() - nRbf, eps); +// Ind += std::pow(Ds,2)*my_pow(delta0,j); + } + } +// Ind *= dist; + + return std::fabs(Ind); +} + +int Multiquadratics_2D_WENO::getDegPoly(int n){ + int degPoly_; + if (n < 5){ + degPoly_ = 1; + }else{ + degPoly_ = (int)floor(-1.5+sqrt(1+8*(n))/2); + } + return degPoly_; +} + diff --git a/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.hpp b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.hpp new file mode 100644 index 0000000..420cdaf --- /dev/null +++ b/ConSol/ConSol/Reconstr/InterpolationBase/Multiquadratics_2D_WENO.hpp @@ -0,0 +1,36 @@ +// +// Created by Fabian Moenkeberg on 2019-12-16. +// + +#ifndef CONSOL_MULTIQUADRATICS_2D_WENO_HPP +#define CONSOL_MULTIQUADRATICS_2D_WENO_HPP + +#include "Multiquadratics_2D.hpp" + +class Multiquadratics_2D_WENO : public Multiquadratics_2D{ +public: + Multiquadratics_2D_WENO(){ + order = 1; + degPoly = 0; + shapeParam = 0.1; + dim = 2; + }; + + Multiquadratics_2D_WENO(int order, double eps); + + Multiquadratics_2D_WENO(int order, int degPoly, double eps); + + ~ Multiquadratics_2D_WENO(){}; + + std::vector evaluate(std::vector> evPoints, MeshBase *mesh, std::vector stencil, std::vector values); + +protected: + std::vector evaluateDirect(std::vector> evPoints, std::vector>> xB, std::vector values, double eps, double rad); + + double ptwiseSmInd(std::vector> xBi, std::vector>> xB, + std::vector values, double eps); + + int getDegPoly(int n); +}; + +#endif //CONSOL_MULTIQUADRATICS_2D_WENO_HPP diff --git a/ConSol/ConSol/Reconstr/MRWENO.cpp b/ConSol/ConSol/Reconstr/MRWENO.cpp new file mode 100644 index 0000000..d87b7cc --- /dev/null +++ b/ConSol/ConSol/Reconstr/MRWENO.cpp @@ -0,0 +1,221 @@ +// +// Created by Fabian Moenkeberg on 2019-12-16. +// + +#include "MRWENO.hpp" + +// +// Created by Fabian Moenkeberg on 30.07.18. +// + +#include "ENO_2D.hpp" +#include +#include "../Utility/MatrixVectorProd.hpp" +//#include "omp.h" +#include +#include +#include "../Utility/ddGeneralizedNewton.hpp" +#include "../Utility/signPropFulfilled.hpp" +#include "../Utility/isWellDef.hpp" +#include "Utility/gauss-legendre-quadrature/gauss_legendre.hpp" +#include "H5Cpp.h" + +void MRWENO::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Configuration config){ + this->model = model; + this->limiter = limiter; + this->interpBase = interpBase; + this->order = config.getOrder(); + this->nThreads = config.getNthreads(); + int n = order; + std::vector w0(n/2 + (n%2)); + std::vector abs(n/2 + (n%2)); + Qlambdas.resize(n, std::vector(2)); + QWeights.resize(n); + double *w = new double[n]; + double *x = new double[n]; + + legendre_set(n,x,w); + if (order < 2){ + Qlambdas[0][0] = 0.5; Qlambdas[0][1]= 0.5; + QWeights[0] = 1; + }else { + for (int i = 0; i < n; i++) { + Qlambdas[i][0] = (1 + x[i]) / 2; + Qlambdas[i][1] = (1 - x[i]) / 2; + QWeights[i] = w[i] / 2; + } + } + + if (nThreads == 0){ + nThreads = omp_get_max_threads(); + } + + std::cout << " Reconstruction: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n"; + + switch (config.getDegPoly()+1) + { + case 2: + nStencil = 5; + maxLevelNeighbour = 3; +// nStencil = 7; +// maxLevelNeighbour = 2; + break; + case 3: + nStencil = 14; + maxLevelNeighbour = 6; +// nStencil = 14; +// maxLevelNeighbour = 3; + break; + case 4: + nStencil = 23; + maxLevelNeighbour = 8; + break; + default: + nStencil = config.getOrder()*(config.getOrder()+1)/2 + 1; + maxLevelNeighbour = ceil(double(nStencil)/2); + break; + } + + delete [] w; + delete [] x; +} + + + +void MRWENO::reconstruct(std::vector>>> &UReconstr, std::vector> &U, MeshBase *mesh, double t){ + int dimSyst = mesh->getNelem(); + int nCells = mesh->getNcells(); + int nQuads = QWeights.size(); + int maxLevel = 3; + std::vector> saveStencil(nCells,std::vector(nStencil)); + + omp_set_num_threads(nThreads); +#pragma omp parallel for + for (int i = 0; i < nCells; i++) { + + int indInternal = mesh->getInternal(i); + std::vector c2edg = mesh->getC2E(indInternal); + double nNodes= c2edg.size(); + std::vector> evPoints(nNodes * nQuads, std::vector(2, 0)); + std::vector> edg_i(2,std::vector(2)); + + for (int l = 0; l < nNodes; l++) { + for (int j = 0; j < nQuads; j++) { + edg_i = mesh->edge2Vertex(c2edg[l]); + evPoints[nNodes * j + l][0] = Qlambdas[j][0] * edg_i[0][0] + Qlambdas[j][1] * edg_i[0][1]; // x-Coordinate + evPoints[nNodes * j + l][1] = Qlambdas[j][0] * edg_i[1][0] + Qlambdas[j][1] * edg_i[1][1]; // y-Coordinate + } + } + + std::vector> uRec(dimSyst, std::vector(nNodes * nQuads)); + for (int k = 0; k < dimSyst; k++) { + int smInd; + for (int level = 0; level < maxLevel; level++) { + std::vector stencil = this->getStencil(indInternal, mesh, &U, k); +// std::vector stencil{1,2,3,4,5}; + int sizeStencil = (int) stencil.size(); + std::vector locValues(sizeStencil); + for (int ii = 0; ii < sizeStencil; ii++) { + locValues[ii] = U[stencil[ii]][k]; + } + std::vector intValues = interpBase->evaluate(evPoints, mesh, stencil, locValues); + uRec[k] = intValues; +// for (int jj = 0; jj < intValues.size(); jj++){ +// uRec[k][jj] = U[indInternal][k]; +// } + } + } + for (int j = 0; j < nNodes; j++) { + if (mesh->getRight(c2edg[j]) == indInternal) { + for (int l = 0; l < nQuads; l++) { + for (int k = 0; k < dimSyst; k++) { + UReconstr[0][c2edg[j]][l][k] = uRec[k][nNodes * l + j]; + } + } + + } else { + for (int l = 0; l < nQuads; l++) { + for (int k = 0; k < dimSyst; k++) { + UReconstr[1][c2edg[j]][l][k] = uRec[k][nNodes * l + j]; + } + } + } + } + } + + mesh->updateEdges(UReconstr, t); +} + +std::vector MRWENO::getStencil(int cellIndex, MeshBase *mesh, std::vector> *U, int k){ + std::vector ret(1); + std::vector ret_tmp(1); + ret[0] = cellIndex; + ret_tmp[0] = cellIndex; + std::vector ids_possible = mesh->getNeightbours(cellIndex); + int nNodes = ids_possible.size(); + for (int kk = 2; kk >= 0; kk--) { + if (ids_possible[kk] == cellIndex) { + ids_possible.erase(ids_possible.begin() + kk); + } + } + std::vector num_possible(ids_possible.size(),1); + std::vector u1(1); + + int maxLevelNeighbour0 = 3; + if (ids_possible.size() < nNodes){ + maxLevelNeighbour0 = std::min(1,maxLevelNeighbour0); + } + u1[0] = U->at(cellIndex)[k]; + +// double smInd0 = interpBase->smoothnessInd(cellIndex, ret_tmp, u1, mesh)/(u1.size()); + double eps0 = 1e-10; + for (int i = 1; i < maxLevelNeighbour0; i++){ + int n = ids_possible.size(); + if (i == 1){ + n = std::min(n, 2); + }else if ( i == 2){ + n = std::min(n, 3); + } + + int min_ind = 0; + + ret_tmp.push_back(ids_possible[0]); + u1.push_back(U->at(ids_possible[0])[k]); + for (int j = 1; j < n; j++){ + ret_tmp.push_back(ids_possible[j]); + u1.push_back(U->at(ids_possible[j])[k]); + } + + for (int j = 0; j < n; j++) { + ret.push_back(ids_possible[0]); + + std::vector newNeighbours = mesh->getNeightbours(ids_possible[j]); + ids_possible.erase(ids_possible.begin() + 0); + + for (int j = 0; j < 3; j++) { + bool addNew = true; + for (int k = 0; k < n - 1; k++) { + if (ids_possible[k] == newNeighbours[j]) { + addNew = false; + } + } + for (int k = 0; k < i + 1; k++) { + if (ret[k] == newNeighbours[j]) { + addNew = false; + } + } + + if (addNew) { + ids_possible.push_back(newNeighbours[j]); + num_possible.push_back(num_possible[min_ind] + 1); + } + } + + num_possible.erase(num_possible.begin() + 0); + } + } + +// std::cout << ret.size() << std::endl; + return ret; +} + diff --git a/ConSol/ConSol/Reconstr/MRWENO.hpp b/ConSol/ConSol/Reconstr/MRWENO.hpp new file mode 100644 index 0000000..579e136 --- /dev/null +++ b/ConSol/ConSol/Reconstr/MRWENO.hpp @@ -0,0 +1,44 @@ +// +// Created by Fabian Moenkeberg on 2019-12-16. +// + +#ifndef CONSOL_MRWENO_HPP +#define CONSOL_MRWENO_HPP + +#include +#include "Reconstruct.hpp" +#include "InterpolationBase/InterpBase.hpp" +#include "Limiter/LimiterBase.hpp" +#include "Utility/Quadrature1D/quadrule.hpp" + +class MRWENO : public Reconstruct{ + int order; + bool reconstrOnCharact = 0; + +protected: + int nStencil; + int maxLevelNeighbour; + +public: + MRWENO(){ + ngc = 1; + name = ENO; + }; + + virtual ~MRWENO(){}; + + void initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Configuration config); + + void reconstruct(std::vector>>> &UReconstr, std::vector> &U, MeshBase *mesh, double t); + +// void exportStencils(std::vector> &U, MeshBase *mesh, double t, std::string output_path); + + +protected: + virtual std::vector getStencil(int cellIndex, MeshBase *mesh, std::vector> *U, int k); + + + void exportVec(std::vector> vec, MeshBase *mesh, std::vector> U, std::string output_path); +}; + +#endif //CONSOL_MRWENO_HPP