diff --git a/ConSol/ConSol/Base/Setup.cpp b/ConSol/ConSol/Base/Setup.cpp index b86755e..fd9a282 100644 --- a/ConSol/ConSol/Base/Setup.cpp +++ b/ConSol/ConSol/Base/Setup.cpp @@ -1,491 +1,495 @@ // // 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: if ((config.getModel()==LINEARADVECTION)||(config.getModel()==BURGERS)||(config.getModel()==KPP)||(config.getModel()==KPPX)||(config.getModel()==KPPY)){ assert(false && "NOt implemented!"); // limiter = ExtremaPresLimiter_2D_TRI(); }else if ((config.getModel()==EULER)){ limiter = new ExtremaPresLimiter_2D_TRI_Euler(config.getOrder(), 1.4); } 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 if (config.getMeshType() == HYBRID){ reconstr = new class HybridRBFWENO(); }else if (config.getStencilChoice() == SYMM) { reconstr = new class CWENO_2D(); }else if (config.getStencilChoice() == SYMMb){ reconstr = new class MRWENO_b(); }else if (config.getStencilChoice() == SYMMc){ reconstr = new class MRWENO_c(); } } break; case NORECONSTR: if (dim ==1){ reconstr = new NoReconstr(); }else if (dim == 2) { - reconstr = new NoReconstr_2D(); + if (config.getMeshType() == HYBRID){ + reconstr = new NoReconstrHybrid_2D(); + }else { + 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_WENO(config.getOrder() - 1, config.getDegPoly(), config.getShapeParam()); // interpBase = new Multiquadratics_2D_WENO_b(config.getOrder() - 1, config.getDegPoly(), // config.getShapeParam()); break; case SYMMb: interpBase = new Multiquadratics_2D_WENO_b(config.getOrder() - 1, config.getDegPoly(), config.getShapeParam()); break; case SYMMc: interpBase = new Multiquadratics_2D_WENO_c(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 NOZZLE: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new NozzleBC(); break; case SCRAMJET: if (dim==1){ assert(false && "Not Implemented, yet!");} else if (config.getInitialCond() == InitNOZZLE){ bc = new ScramJetBC(0.5); } else bc = new ScramJetBC(); break; case SCRAMJET10: if (dim==1){ assert(false && "Not Implemented, yet!");} else bc = new ScramJetBC(10.0); 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; case AXISYMMEULER: source = new AxisymmEuler_std_noSwirl(); break; case AXISYMMEULERRU: source = new AxisymmEuler_rU_noSwirl(); 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 if (config.getBc() == SCRAMJET10){ initialCond = new initFAC(0, 10.0); }else { initialCond = new initFAC(); } break; case InitAIRFOIL: initialCond = new AirFoil(); break; case InitNOZZLE: initialCond = new initNozzle(); break; case InitAXISYMMNOZZLE: initialCond = new initAxiSymmNozzle(); break; case InitAXISYMMNOZZLERU: initialCond = new initAxiSymmNozzle_rU(); 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); std::cout << "Mesh initialized" << std::endl; interpBase->initialize(config); std::cout << "Interpolation Base initialized" << std::endl; if (source->getIsUsed()){ source->initialize(config, mesh); std::cout << "Source term initialized" << std::endl; } reconstr->initialize(model, interpBase, limiter, mesh, source, config); std::cout << "Reconstruction initialized" << std::endl; flux->initialize(mesh, source, model, reconstr, config); std::cout << "Flux initialized" << std::endl; limiter->initialize(interpBase, model); std::cout << "Limiter initialized" << std::endl; timeInt->initialize(config, mesh); std::cout << "Time integration initialized" << std::endl; this->soln = Solution(config, mesh, initialCond); std::cout << "Solution initialized" << std::endl; 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 e5cf1c2..5affb52 100644 --- a/ConSol/ConSol/Base/Setup.hpp +++ b/ConSol/ConSol/Base/Setup.hpp @@ -1,215 +1,216 @@ // // 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 "../Source/AxisymmEuler_std_noSwirl.hpp" #include "../Source/AxisymmEuler_rU_noSwirl.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_2D_TRI.hpp" #include "../Reconstr/Limiter/ExtremaPresLimiter_2D_TRI_Euler.hpp" #include "../Reconstr/NoReconstr.hpp" #include "../Reconstr/NoReconstr_2D.hpp" +#include "../Reconstr/NoReconstrHybrid_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 "../Utility/InitialConditions/2D/InitNozzle.hpp" #include "../Utility/InitialConditions/2D/InitAxiSymmNozzle.hpp" #include "../Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.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/Multiquadratics_2D_WENO_b.hpp" #include "../Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.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/MRWENO_b.hpp" #include "../Reconstr/MRWENO_c.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" #include "../Mesh/BC/NozzleBC.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 34f7de5..3c81bd2 100644 --- a/ConSol/ConSol/CMakeLists.txt +++ b/ConSol/ConSol/CMakeLists.txt @@ -1,314 +1,314 @@ 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/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 Reconstr/MRWENO.cpp Reconstr/MRWENO.hpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO.cpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO.hpp Utility/InitialConditions/2D/InitNozzle.cpp Utility/InitialConditions/2D/InitNozzle.hpp Mesh/BC/flux_inflow.cpp Mesh/BC/flux_inflow.hpp Mesh/BC/flux_outflow.cpp - Mesh/BC/flux_outflow.hpp Mesh/BC/NozzleBC.cpp Mesh/BC/NozzleBC.hpp Model/EulerAxiSymm.cpp Model/EulerAxiSymm.hpp Source/AxisymmEuler.cpp Source/AxisymmEuler.hpp Reconstr/MRWENO_b.cpp Reconstr/MRWENO_b.hpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_b.cpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_b.hpp Reconstr/MRWENO_c.cpp Reconstr/MRWENO_c.hpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.cpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.hpp Source/AxisymmEuler_std_noSwirl.cpp Source/AxisymmEuler_std_noSwirl.hpp Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp Utility/InitialConditions/2D/InitAxiSymmNozzle.hpp Source/AxisymmEuler_rU_noSwirl.cpp Source/AxisymmEuler_rU_noSwirl.hpp Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.cpp Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.hpp Reconstr/InterpolationBase/Multiquadratics_2D_CWENO.cpp Reconstr/InterpolationBase/Multiquadratics_2D_CWENO.hpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI.cpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI.hpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI_Euler.cpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI_Euler.hpp Reconstr/Limiter/ExtremaPresLimiter_2D_QUAD.cpp Reconstr/Limiter/ExtremaPresLimiter_2D_QUAD.hpp) + Mesh/BC/flux_outflow.hpp Mesh/BC/NozzleBC.cpp Mesh/BC/NozzleBC.hpp Model/EulerAxiSymm.cpp Model/EulerAxiSymm.hpp Source/AxisymmEuler.cpp Source/AxisymmEuler.hpp Reconstr/MRWENO_b.cpp Reconstr/MRWENO_b.hpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_b.cpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_b.hpp Reconstr/MRWENO_c.cpp Reconstr/MRWENO_c.hpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.cpp Reconstr/InterpolationBase/Multiquadratics_2D_WENO_c.hpp Source/AxisymmEuler_std_noSwirl.cpp Source/AxisymmEuler_std_noSwirl.hpp Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp Utility/InitialConditions/2D/InitAxiSymmNozzle.hpp Source/AxisymmEuler_rU_noSwirl.cpp Source/AxisymmEuler_rU_noSwirl.hpp Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.cpp Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.hpp Reconstr/InterpolationBase/Multiquadratics_2D_CWENO.cpp Reconstr/InterpolationBase/Multiquadratics_2D_CWENO.hpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI.cpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI.hpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI_Euler.cpp Reconstr/Limiter/ExtremaPresLimiter_2D_TRI_Euler.hpp Reconstr/Limiter/ExtremaPresLimiter_2D_QUAD.cpp Reconstr/Limiter/ExtremaPresLimiter_2D_QUAD.hpp Reconstr/NoReconstrHybrid_2D.cpp Reconstr/NoReconstrHybrid_2D.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/NoReconstrHybrid_2D.cpp b/ConSol/ConSol/Reconstr/NoReconstrHybrid_2D.cpp new file mode 100644 index 0000000..e4e873f --- /dev/null +++ b/ConSol/ConSol/Reconstr/NoReconstrHybrid_2D.cpp @@ -0,0 +1,144 @@ +// +// Created by Fabian Moenkeberg on 2020-05-26. +// + +#include "NoReconstrHybrid_2D.hpp" + +void NoReconstrHybrid_2D::initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Source* source, Configuration config){ + + this->model = model; + this->source = source; + this->limiter = limiter; + this->interpBase = interpBase; + this->order = ceil((config.getOrder()+1)/2.0)+1; + this->nThreads = config.getNthreads(); + + if (nThreads == 0){ + nThreads = omp_get_max_threads(); + } + + std::cout << " Reconstruction: Running threads: " << nThreads << " out of " << omp_get_max_threads() << "\n"; + + xCoord = mesh->getXCoord(); + ngc = mesh->getNgc(); + + + nCells = mesh->getNcells(); + nQuads = QWeights.size(); + nDomainsQUAD = mesh->getNdomainsQUAD() + mesh->getNdomainsQUADQUAD() + mesh->getNdomainsRQ(); + + nDomainsTRIQUAD = mesh->getNdomainsTRIQUAD(); + nDomainsTRI = mesh->getNdomainsTRI(); + nDomainsRT = mesh->getNdomainsRT(); + submeshCell0 = mesh->getSubmeshCell0(); + submeshEdg0 = mesh->getSubmeshEdg0(); + nX = mesh->getNx(); + nY = mesh->getNy(); +} + +void NoReconstrHybrid_2D::reconstruct(std::vector>>> &UReconstr, std::vector> &U, + MeshBase *mesh, double t, std::vector> &Source){ + + int dimSyst = mesh->getNelem(); + int nCells = mesh->getNcells(); + UReconstr.resize(2,std::vector>>(1,std::vector>(mesh->getNedges(),std::vector(dimSyst)))); + std::vector source_loc; + + std::vector c2e_loc(4); + int nNodes; + for (int ii = 0; ii < nCells; ii++){ + c2e_loc = mesh->getC2E_internal(ii); + nNodes = c2e_loc.size(); + int internal_i = mesh->getInternal(ii); + + for (int j = 0; j < nNodes; j++){ + + if (mesh->getRight(c2e_loc[j])== internal_i){ + UReconstr[0][c2e_loc[j]][0] = U[internal_i]; + }else{ + UReconstr[1][c2e_loc[j]][0] = U[internal_i]; + } + } + if (source->getIsUsed()) { + std::vector> xy = mesh->getXCoordOfCell(internal_i); + source_loc = source->source(U[internal_i], t, std::vector{(xy[0][0] + xy[1][0] + xy[2][0]) / 3, + (xy[0][1] + xy[1][1] + xy[2][1]) / 3}); + + for (int k = 0; k < dimSyst; k++) { + Source[internal_i][k] = source_loc[k]; + } + } + } + mesh->updateEdges(UReconstr, t); + + for (int ii = 0; ii < nDomainsQUAD; ii++){ + omp_set_num_threads(nThreads); +#pragma omp parallel for + for (int i = 0; i < nX[ii]; i++) { + std::vector> uRec(dimSyst, std::vector(2)); + std::vector uloc(2*(order-1) + 1); + std::vector rec(2); + + for (int j = 0; j < nY[ii]-1; j++) { + for (int k = 0; k < dimSyst; k++) { + uRec[k][0] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + uRec[k][1] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + } + + limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]); + for (int k = 0; k < dimSyst; k++) { + UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0]; + UReconstr[1][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][1]; + } + } + int j = nY[ii]-1; + + for (int k = 0; k < dimSyst; k++) { + uRec[k][0] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + uRec[k][1] = U[i+ngc+(j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + } + + limiter->limit_WENO(uRec, 0, mesh, U[i+ngc+(0+j+ngc)*(nX[ii]+2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]); + for (int k = 0; k < dimSyst; k++) { + UReconstr[0][i+(j)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][0]; + UReconstr[0][i+(j+1)*nX[ii] + submeshEdg0[ii + 1]][0][k] = uRec[k][1]; + } + } + omp_set_num_threads(nThreads); +#pragma omp parallel for + for (int j = 0; j < nY[ii]; j++) { + std::vector> uRec(dimSyst, std::vector(2)); + std::vector uloc(2*(order-1) + 1); + std::vector rec(2); + + for (int i = 0; i < nX[ii]-1; i++) { + + for (int k = 0; k < dimSyst; k++) { + uRec[k][0] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + uRec[k][1] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + } + + limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]); + for (int k = 0; k < dimSyst; k++) { + UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0]; + UReconstr[1][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1]; + } + } + + int i = nX[ii] - 1; + for (int k = 0; k < dimSyst; k++) { + uRec[k][0] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + uRec[k][1] = U[i + ngc + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD]][k]; + } + + limiter->limit_WENO(uRec, 0, mesh, U[i + ngc + 0 + (j+ngc)*(nX[ii] + 2*ngc) + submeshCell0[ii + nDomainsTRI + nDomainsTRIQUAD + nDomainsRT]]); + for (int k = 0; k < dimSyst; k++) { + UReconstr[0][i+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][0]; + UReconstr[0][i+1+(j)*(nX[ii] + 1) + nX[ii]*(nY[ii] + 1) + submeshEdg0[ii + 1]][0][k] = uRec[k][1]; + } + } + } + + mesh->updateEdges(UReconstr, t); +} + diff --git a/ConSol/ConSol/Reconstr/NoReconstrHybrid_2D.hpp b/ConSol/ConSol/Reconstr/NoReconstrHybrid_2D.hpp new file mode 100644 index 0000000..f3232fb --- /dev/null +++ b/ConSol/ConSol/Reconstr/NoReconstrHybrid_2D.hpp @@ -0,0 +1,47 @@ +// +// Created by Fabian Moenkeberg on 2020-05-26. +// + +#ifndef CONSOL_NORECONSTRHYBRID_2D_HPP +#define CONSOL_NORECONSTRHYBRID_2D_HPP + + +#include "Reconstruct.hpp" + +class NoReconstrHybrid_2D : public Reconstruct{ + int order; + std::vector> xCoord; + std::vector> Crec; + std::vector nY; + std::vector nX; + std::vector submeshCell0; + std::vector submeshEdg0; + int nCells; + int nQuads; + int nDomainsQUAD; + int nDomainsTRI; + int nDomainsTRIQUAD; + int nDomainsRT; + + int order_scheme; +public: + NoReconstrHybrid_2D(){ + ngc = 1; + name = NORECONSTR; + this->QWeights.resize(1); + this->QWeights[0] = 1.0; + }; + + ~NoReconstrHybrid_2D(){}; + + NoReconstrHybrid_2D(ModelBase *model, Configuration config); + + void initialize(ModelBase *model, InterpBase *interpBase, LimiterBase *limiter, MeshBase *mesh, Source* source, Configuration config); + +// virtual void reconstruct(std::vector>> &UReconstr, std::vector> U, MeshBase *mesh) ; + + virtual void reconstruct(std::vector>>> &UReconstr, std::vector> &U, + MeshBase *mesh, double t, std::vector> &Source); +}; + +#endif //CONSOL_NORECONSTRHYBRID_2D_HPP diff --git a/ConSol/Examples/HybridTest/Airfoil/CalcPressureAtFoil/calcPressureDistr.m b/ConSol/Examples/HybridTest/Airfoil/CalcPressureAtFoil/calcPressureDistr.m index 2cb0052..a6a5ac5 100644 --- a/ConSol/Examples/HybridTest/Airfoil/CalcPressureAtFoil/calcPressureDistr.m +++ b/ConSol/Examples/HybridTest/Airfoil/CalcPressureAtFoil/calcPressureDistr.m @@ -1,61 +1,61 @@ % get pressure profile a=importdata('NACA.txt'); filename = '../Results/Airfoil_otf28_NACA0012_190508_frontd_t150_newBC.h5'; hinfo = hdf5info(filename); dim = hinfo.GroupHierarchy.Groups.Attributes(1).Value; nSyst = hinfo.GroupHierarchy.Groups.Attributes(2).Value; timeElapsed = hinfo.GroupHierarchy.Attributes.Value cells = hdf5read(hinfo.GroupHierarchy.Groups.Datasets(1))+1; internal = hdf5read(hinfo.GroupHierarchy.Groups.Datasets(2))+1; xCoord = hdf5read(hinfo.GroupHierarchy.Groups.Datasets(3)); time = hdf5read(hinfo.GroupHierarchy.Datasets(1)); rho = hdf5read(hinfo.GroupHierarchy.Datasets(1+1)); soln = rho; m1 = hdf5read(hinfo.GroupHierarchy.Datasets(2+1)); m2 = hdf5read(hinfo.GroupHierarchy.Datasets(3+1)); E = hdf5read(hinfo.GroupHierarchy.Datasets(4+1)); p = (1.4 - 1) * (E - 0.5*(m1.*m1 + m2.*m2) ./ rho); M = sqrt((m1.^2+m2.^2)./(1.4*p.*rho)); zLim = [min(M(:,end)),max(M(:,end))]; zLim2 = [min(rho(:,end)),max(rho(:,end))]; area = 1/2*abs((xCoord(cells(2,:),1)-xCoord(cells(1,:),1)).*(xCoord(cells(3,:),2)... -xCoord(cells(1,:),2))-(xCoord(cells(3,:),1)-xCoord(cells(1,:),1)).*(xCoord(cells(2,:),2)-xCoord(cells(1,:),2))); % xM = (xCoord(cells(1,:),:) + xCoord(cells(2,:),:) + xCoord(cells(3,:),:) + xCoord(cells(4,:),:))/4; xM = (xCoord(cells(1,:),:) + xCoord(cells(2,:),:) + xCoord(cells(3,:),:))/3; a.data = (a.data(1:end-1,:)+a.data(2:end,:))/2; N0 = size(a.data,1); pP = zeros(N0,1); IND = zeros(N0,1); for i = 1:N0 [value,ind] = min(sqrt((xM(:,1)-a.data(i,1)).^2 + (xM(:,2)-a.data(i,2)).^2)); - pP(i) = p(ind); + pP(i) = p(ind,end); % if ind == 7690 % ind % end IND(i) = ind; end Cp = 2/1.4/0.85^2*(pP/1-1); -figure(11); +figure(10); plot(a.data(1:N0/2,1),Cp(1:N0/2),'b'); hold on; plot(a.data(N0/2+1:end,1),Cp(N0/2+1:end),'r'); set(gca,'Ydir','reverse') xlabel('x'); ylabel('C_p'); figure(12); plot(a.data(1:N0/2,1),a.data(1:N0/2,2)); hold on;plot(a.data(N0/2+1:end,1),a.data(N0/2+1:end,2)); for i = 1:N0 plot(xM(IND(i),1),xM(IND(i),2),'b*'); end \ No newline at end of file