diff --git a/ConSol/ConSol/Base/Configuration.cpp b/ConSol/ConSol/Base/Configuration.cpp index 2118ff2..aa65f45 100644 --- a/ConSol/ConSol/Base/Configuration.cpp +++ b/ConSol/ConSol/Base/Configuration.cpp @@ -1,896 +1,900 @@ // // Configuration.cpp // ConserSolver // // Created by Fabian Moenkeberg on 04.05.17. // Copyright © 2017 FabianMoenkeberg. All rights reserved. // #include "Configuration.hpp" #include #include Configuration::Configuration(){ }; Configuration::Configuration(char* file) { this->input_file = file; read(); this->SOLN_DIR = OUTPUT_PATH + "/SOLN_FILES"; // // Making solution directories // std::string commandline0 = "rm -rf " + SOLN_DIR; // system(commandline0.c_str()); // std::string commandline1 = "mkdir -p " + SOLN_DIR; // system(commandline1.c_str()); // // glob_filename = SOLN_DIR + "/globals.dat"; // glob_file.open (glob_filename.c_str()); // // L1_err_filename = SOLN_DIR + "/l1_err.dat"; // L2_err_filename = SOLN_DIR + "/l2_err.dat"; // Linf_err_filename = SOLN_DIR + "/linf_err.dat"; // // if(param.exact_available) // { // L1_err_file.open (L1_err_filename.c_str()); // L2_err_file.open (L2_err_filename.c_str()); // Linf_err_file.open (Linf_err_filename.c_str()); // } // // if(param.OOC) // { // print_globals = false; // print_soln = false; // disp_log = false; // L1_err.resize(param.mesh_levels); // L2_err.resize(param.mesh_levels); // Linf_err.resize(param.mesh_levels); // } // else { // print_globals = true; // print_soln = true; // disp_log = true; // } } Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, double cfl){ this->model = model; this->evalMethod = evalMethod; this->bcname = bc; this->order = order; this->tMax = tMax; this->cfl = cfl; this->tInclude = *new std::vector(0); this->reconstr = NORECONSTR; this->interpBase = POLYNOMIALS; } Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, std::vector tInclude, double cfl){ this->model = model; this->evalMethod = evalMethod; this->bcname = bc; this->order = order; this->tMax = tMax; this->cfl = cfl; this->tInclude = tInclude; this->reconstr = NORECONSTR; this->interpBase = POLYNOMIALS; } Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, double cfl){ this->model = model; this->evalMethod = evalMethod; this->bcname = bc; this->order = order; this->tMax = tMax; this->cfl = cfl; this->tInclude = *new std::vector(0); this->reconstr = reconstr; this->interpBase = interpBase; } Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, std::vector tInclude, double cfl){ this->model = model; this->evalMethod = evalMethod; this->bcname = bc; this->order = order; this->tMax = tMax; this->cfl = cfl; this->tInclude = tInclude; this->reconstr = reconstr; this->interpBase = interpBase; } Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, double cfl){ this->model = model; this->evalMethod = evalMethod; this->bcname = bc; this->order = order; this->tMax = tMax; this->cfl = cfl; this->tInclude = *new std::vector(0); this->reconstr = reconstr; this->interpBase = interpBase; this->reconstrOnCharact = reconstrOnCharact; } Configuration::Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, std::vector tInclude, double cfl){ this->model = model; this->evalMethod = evalMethod; this->bcname = bc; this->order = order; this->tMax = tMax; this->cfl = cfl; this->tInclude = tInclude; this->reconstr = reconstr; this->interpBase = interpBase; this->reconstrOnCharact = reconstrOnCharact; } ModelName Configuration::getModel(){return model;}; TimeIntegrationName Configuration::getTimeIntegration(){return time_scheme;}; int Configuration::getOrder(){return order;}; evaluationMethod Configuration::getEvalMethod(){return evalMethod;}; bcName Configuration::getBc(){return bcname;}; std::vector Configuration::getTInclude(){return tInclude;}; limiterName Configuration::getLimiter(){return limiter_name; }; FluxName Configuration::getFlux(){return fluxScheme; }; double Configuration::getTmax(){return tMax;}; double Configuration::getCfl(){return cfl;}; int Configuration::getMaxNumWrite(){return maxNumWrite;}; ReconstrName Configuration::getReconstr(){return reconstr;}; sourceName Configuration::getSource() { return source;}; initalCondName Configuration::getInitialCond() {return initialCond;}; InterpBaseName Configuration::getInterpBase(){return interpBase;}; bool Configuration::getReconstrOnCharact(){return reconstrOnCharact;}; int Configuration::getDim(){return dim;}; double Configuration::getXmax(){return xmax;}; double Configuration::getXmin(){return xmin;}; double Configuration::getYmax(){return ymax;}; double Configuration::getYmin(){return ymin;}; double Configuration::getShapeParam(){ return shapeParam;}; std::vector Configuration::getNx() { return Nx_cell;}; std::vector Configuration::getNy() { return Ny_cell;}; std::string Configuration::getGridFile() { return grid_file;}; std::string Configuration::getOutputPath() { return OUTPUT_PATH;} std::string Configuration::getImportPath() { return IMPORT_PATH;} bool Configuration::getImportSol() {return import_sol;}; bool Configuration::getLoadGrid() { return load_grid;}; int Configuration::getNthreads() {return nThreads;}; stencilChoice Configuration::getStencilChoice() { return stencilchoice;}; std::map Configuration::getConstants() {return constants;}; int Configuration::getDegPoly() {return degPoly;}; MeshType Configuration::getMeshType() {return meshType;}; //------------------------------------------------------------------------------ // Read parameters from file //------------------------------------------------------------------------------ void Configuration::read () { std::cout << "Reading input file " << input_file << std::endl; Reader fin(input_file); read_grid (fin); read_numeric (fin); read_model(fin); read_constants (fin); read_initial_condition (fin); read_exact (fin); read_boundary_condition (fin); read_output (fin); read_import(fin); } //------------------------------------------------------------------------------ // Read grid section //------------------------------------------------------------------------------ void Configuration::read_grid (Reader &fin) { std::cout << " Reading grid section\n"; std::string input; int N; fin.begin_section ("grid"); fin.entry ("load_grid"); fin >> input; if (input == "true") { load_grid = true; } else{ load_grid = false; } fin.entry ("dim"); fin >> dim; assert(dim == 1 || dim == 2 ); fin.entry ("x_min"); fin >> xmin; fin.entry ("x_max"); fin >> xmax; fin.begin_section ("Nx_cell"); while(!fin.eos()) { fin >> N; assert(N > 0); Nx_cell.push_back(N); } if (dim == 2) { fin.entry("y_min"); fin >> ymin; fin.entry("y_max"); fin >> ymax; }else{ ymin = 0; ymax = 1; } fin.begin_section ("Ny_cell"); while(!fin.eos()) { fin >> N; assert(N > 0); Ny_cell.push_back(N); } assert(xmin> grid_file; fin.entry ("type"); fin >> input; if (input == "triangular") { meshType = TRI; } else if (input == "hybrid") { meshType = HYBRID; }else if (input == "quadrilateral"){ meshType = QUAD; }else if (input == "mixed"){ meshType = MIXED; } else { std::cout << " Error: unknown meshType " << input << std::endl; exit(0); } fin.end_section (); } //------------------------------------------------------------------------------ // Read numeric section //------------------------------------------------------------------------------ void Configuration::read_numeric (Reader &fin) { std::cout << " Reading numeric section\n"; std::string input; fin.begin_section("numeric"); fin.entry("nthreads"); fin >> nThreads; fin.entry("time_scheme"); fin >> input; if (input == "SSPRK3") { time_scheme = SSPRK3; } else if (input == "SSPRK5") { time_scheme = SSPRK5; }else if (input == "SSPRK2"){ time_scheme = SSPRK2; }else if (input == "FE") { time_scheme = FE; } else { std::cout << " Error: unknown time_scheme " << input << std::endl; exit(0); } fin.entry("cfl"); fin >> cfl; assert (cfl > 0.0); fin.entry("max_iter"); fin >> max_iter; assert (max_iter > 0); fin.entry("final_time"); fin >> tMax; assert (tMax >= 0.0); fin.entry("min_residue"); fin >> min_residue; assert (min_residue > 0.0); fin.entry("reconstruct"); fin >> input; if (input == "noReconstruct") { reconstr = NORECONSTR; } else if (input == "eno") { reconstr = ENO; } else if (input == "weno") { reconstr = WENO; } else { std::cout << " Error: unknown reconstruction method " << input << std::endl; exit(0); } fin.entry("stencilChoice"); fin >> input; if (input == "std") { stencilchoice = STD; } else if (input == "otf") { stencilchoice = OTF; } else if (input == "symm") { stencilchoice = SYMM; } else if (input == "symmb") { stencilchoice = SYMMb; } else if (input == "symmc") { stencilchoice = SYMMc; } else { std::cout << " Error: unknown stencil choice " << input << std::endl; exit(0); } fin.entry("order"); fin >> order; fin.entry("degPoly"); fin >> input; if (input == "x"){ degPoly = order - 2; } else { degPoly = std::stoi(input); } fin.entry("interpolationBase"); fin >> input; if (input == "Multiquadratics") { interpBase = MULTIQUADS; } else if (input == "Polynomials") { interpBase = POLYNOMIALS; } else if (input == "Polyharmonics") { interpBase = POLYHARMONICS; } else if (input == "Gaussians") { interpBase = GAUSSIANS; } else { std::cout << " Error: unknown interpolation base " << input << std::endl; exit(0); } fin.entry("shapeParameter"); fin >> shapeParam; fin.entry("reconstrOnCharact"); fin >> input; if (input == "true") { reconstrOnCharact = true; } else{ reconstrOnCharact = false; } fin.entry("evaluationMethod"); fin >> input; if (input == "meanvalue") { evalMethod = MEANVALUE; } else{ evalMethod = PTWISE; } fin.entry("fluxScheme"); fin >> input; if (input == "LaxFried") { fluxScheme = LAXFRIEDRICHS; } else { std::cout << " Error: unknown flux scheme " << input << std::endl; exit(0); } fin.entry("limiter"); fin >> input; if (input == "noLimiter") { limiter_name = NOLIMITER; } else if (input == "extremaPreserving") { limiter_name = EXTREMAPRES; } else { std::cout << " Error: unknown limiter " << input << std::endl; exit(0); } fin.begin_section("tInclude"); double value; this->tInclude = std::vector(0); while (!fin.eos()) { fin >> value; tInclude.push_back(value); } fin.end_section (); } //------------------------------------------------------------------------------ // Read material section //------------------------------------------------------------------------------ void Configuration::read_model (Reader &fin) { std::cout << " Reading model section\n"; std::string input; fin.begin_section ("model"); fin.entry("model"); fin >> input; if (input == "linAdv") { model = LINEARADVECTION; } else if (input == "burgers") { model = BURGERS; } else if (input == "SW") { model = SW; } else if (input == "Euler") { model = EULER; } else if (input == "BuckleyLeverett") { model = BUCKLEYLEVERETT; }else if (input == "KPP"){ model = KPP; }else if (input == "KPPY"){ model = KPPY; }else if (input == "KPPX"){ model = KPPX; } else { std::cout << " Error: unknown model " << input << std::endl; exit(0); } fin.entry("source"); fin >> input; if (input == "noSource") { source = NOSOURCE; }else if (input == "axisymmEuler") { source = AXISYMMEULER; + }else if (input == "axisymmEulerRU") { + source = AXISYMMEULERRU; }else{ std::cout << " Error: unknown source " << input << std::endl; exit(0); } fin.end_section (); - if (source == AXISYMMEULER && model!=EULER){ + if ((source == AXISYMMEULER|| source == AXISYMMEULERRU) && model!=EULER){ assert(false && "Configuration: Axisymmetric euler source without euler model!"); } } //------------------------------------------------------------------------------ // Read constants section //------------------------------------------------------------------------------ void Configuration::read_constants (Reader &fin) { std::cout << " Reading constants section\n"; std::string input; double value; fin.begin_section ("constants"); while (!fin.eos()) { fin >> input; fin >> value; std::cout << std::setw(16) << input << std::setw(16) << value << std::endl; constants.insert ( std::pair(input, value) ); } } //------------------------------------------------------------------------------ // Read initial condition //------------------------------------------------------------------------------ void Configuration::read_initial_condition (Reader &fin) { std::cout << " Reading initial condition section\n"; std::string input; fin.begin_section ("initial_condition"); switch (model){ case LINEARADVECTION: fin.entry("y"); fin >> input; if (input == "sin(2*pi*x)"){ initialCond = SIN; }else if (input == "bubble") { initialCond = BUBBLE; } else if ( input == "coscos") { initialCond = COSCOS; } else if ( input == "sinsin") { initialCond = SINSIN; } else if (input == "shockX"){ initialCond = SHOCKX; } else if (input == "cube"){ initialCond = InitCUBE; }else { std::cout << " Error:: unknown initial Condition: " << input << std::endl; exit (0); } break; case EULER: fin.entry("y"); fin >> input; if (input == "isentropicVortex"){ initialCond = InitISENTROPICVORTEX; }else if (input == "eulerBlast"){ initialCond = InitEULERBLAST; }else if (input == "vortexShock"){ initialCond = InitVORTEXSHOCK; }else if (input == "RP"){ initialCond = RP; }else if (input == "initDoubleMachRef"){ initialCond = InitDOUBLEMACHREF; }else if (input == "initFAC"){ initialCond = InitFAC; }else if (input == "initAirFoil"){ initialCond = InitAIRFOIL; }else if (input == "initNozzle"){ initialCond = InitNOZZLE; }else if (input == "initAxisymmNozzle"){ initialCond = InitAXISYMMNOZZLE; + }else if (input == "initAxisymmNozzleRU"){ + initialCond = InitAXISYMMNOZZLERU; } else { std::cout << " Error:: unknown initial Condition: " << input << std::endl; exit (0); } break; case BURGERS: fin.entry("y"); fin >> input; if (input == "sin(2*pi*x)"){ initialCond = SIN; }else if (input == "bubble") { initialCond = BUBBLE;} else if (input == "burgers") { initialCond = InitBURGERS;} else{ std::cout << " Error:: unknown initial Condition: " << input << std::endl; exit (0); } break; case SW: std::cout << " Error:: Not implemented, yet: " << input << std::endl; exit (0); break; case BUCKLEYLEVERETT: std::cout << " Error:: Not implemented, yet: " << input << std::endl; exit (0); break; case KPP: fin.entry("y"); fin >> input; if (input == "KPP") { initialCond = InitKPP; }else { std::cout << " Error:: unknown initial Condition: " << input << std::endl; exit (0); } break; case KPPY: fin.entry("y"); fin >> input; if (input == "KPP") { initialCond = InitKPP; }else if (input == "KPPY") { initialCond = InitKPPY; }else { std::cout << " Error:: unknown initial Condition: " << input << std::endl; exit (0); } break; case KPPX: fin.entry("y"); fin >> input; if (input == "KPP") { initialCond = InitKPP; }else { std::cout << " Error:: unknown initial Condition: " << input << std::endl; exit (0); } break; } // // // fin.entry("density"); // fin.getline (input); // initial_condition.add ("density", input); // // fin.entry ("xvelocity"); // fin.getline (input); // initial_condition.add ("xvelocity", input); // // fin.entry ("yvelocity"); // fin.getline (input); // initial_condition.add ("yvelocity", input); // // fin.entry ("zvelocity"); // fin.getline (input); // initial_condition.add ("zvelocity", input); // // fin.entry ("pressure"); // fin.getline (input); // initial_condition.add ("pressure", input); // // fin.end_section (); } //------------------------------------------------------------------------------ // Read exact solution //------------------------------------------------------------------------------ void Configuration::read_exact (Reader &fin) { std::cout << " Reading exact solution section\n"; std::string input; fin.begin_section ("exact_soln"); fin.entry("available"); fin >> input; if(input == "yes") exact_available = true; else if(input == "no") exact_available = false; else { std::cout << " Error:: unknown exact solution availability option: " << input << std::endl; exit (0); } if(exact_available) { // fin.entry("density"); // fin.getline (input); // exact_soln.add ("density", input); // // fin.entry ("xvelocity"); // fin.getline (input); // exact_soln.add ("xvelocity", input); // // fin.entry ("yvelocity"); // fin.getline (input); // exact_soln.add ("yvelocity", input); // // fin.entry ("zvelocity"); // fin.getline (input); // exact_soln.add ("zvelocity", input); // // fin.entry ("pressure"); // fin.getline (input); // exact_soln.add ("pressure", input); // // fin.end_section (); } else { while(!fin.eos()) fin >> input; } } //------------------------------------------------------------------------------ // Read boundary conditions //------------------------------------------------------------------------------ void Configuration::read_boundary_condition (Reader &fin) { std::cout << " Reading boundary condition section\n"; std::string input; fin.begin_section ("boundary_condition"); fin.entry("bc"); fin >> input; if(input == "mixedBC") bcname = MIXEDBC; else if(input == "periodic") bcname = PERIODIC; else if(input == "absorb") bcname = ABSORBINGBC; else if(input == "slipWall") bcname = SLIPWALL; else if(input == "doubleMachRefBC") bcname = DOUBLEMACHREF; else if(input == "FAC") bcname = FAC; else if(input == "OpenFlow") bcname = OPENFLOW; else if(input == "ScramJet") bcname = SCRAMJET; else if(input == "Nozzle") bcname = NOZZLE; else { std::cout << " Error:: unknown boundary condition: " << input << std::endl; exit (0); } std::string boundary_names[] = {"xleft","xright","yleft","yright","zleft","zright"}; for(int i=0; i<2*dim;++i) { fin.begin_section (boundary_names[i]); fin.entry("type"); fin >> input; if(input=="periodic") { bcType[i] = PERIODIC; } else if(input=="dirichlet") { bcType[i] = DIRICHLET; } else if(input=="donothing") { bcType[i] = NEUMANN; } else if(input=="slipWall") { bcType[i] = SLIPWALL; } else if(input == "doubleMachRefBC") bcType[i] = DOUBLEMACHREF; else if(input == "absorb") bcType[i] = ABSORBINGBC; else if(input == "FAC") bcType[i] = FAC; else if(input == "OpenFlow") bcType[i] = OPENFLOW; else if(input == "ScramJet") bcType[i] = SCRAMJET; else if(input == "Nozzle") bcType[i] = NOZZLE; else { std::cout << " Error: unknown boundary type " << input << std::endl; exit (0); } if(input == "dirichlet") { std::cout << " Error: Dirichlet type BC not implemented yet "<< std::endl; exit (0); } fin.end_section(); } fin.end_section (); } //------------------------------------------------------------------------------ // Read output section //------------------------------------------------------------------------------ void Configuration::read_output (Reader &fin) { std::cout << " Reading output section\n"; std::string input; fin.begin_section ("output"); // fin.entry("OOC_study"); // fin >> input; // if(input == "yes") // OOC = true; // else if(input == "no") // OOC = false; // else // { // std::cout << " Error:: unknown OOC_study option: " << input << std::endl; // exit (0); // } // if(OOC && (Nx_cell.size() <2 || Ny_cell.size() < 2)) // { // std::cout << " Error:: Cannot perform OOC study without at least two mesh levels"<< std::endl; // exit (0); // } // // if(OOC && !exact_available) // { // std::cout << " Error:: Cannot perform OOC study as exact solution not avaible"<> t_stamps; // assert (t_stamps >= 1); fin.entry ("output_path"); fin >> OUTPUT_PATH; fin.end_section (); } //------------------------------------------------------------------------------ // Read Import section //------------------------------------------------------------------------------ void Configuration::read_import (Reader &fin) { std::cout << " Reading import section\n"; std::string input; fin.begin_section ("import"); fin.entry ("import"); fin >> input; if (input == "true") { import_sol = true; } else{ import_sol = false; } fin.entry ("import_path"); fin >> IMPORT_PATH; fin.end_section (); } \ No newline at end of file diff --git a/ConSol/ConSol/Base/Configuration.hpp b/ConSol/ConSol/Base/Configuration.hpp index 74dc59c..8ae13ef 100644 --- a/ConSol/ConSol/Base/Configuration.hpp +++ b/ConSol/ConSol/Base/Configuration.hpp @@ -1,189 +1,189 @@ // // Configuration.hpp // ConserSolver // // Created by Fabian Moenkeberg on 04.05.17. // Copyright © 2017 FabianMoenkeberg. All rights reserved. // #ifndef Configuration_hpp #define Configuration_hpp //class ModelBase; #include #include #include #include #include "reader.h" #include enum evaluationMethod{PTWISE, MEANVALUE}; enum bcName{INFLOW, OUTFLOW, SLIPWALL, EXACT, PERIODIC, DIRICHLET, NEUMANN, MIXEDBC, MIXEDBCCOMPONENT, SPECIAL, DOUBLEMACHREF , Q2I, T2I, ABSORBINGBC, Q2Q, FAC, OPENFLOW, SCRAMJET, NOZZLE}; enum ModelName{EULER, BURGERS,SW,BUCKLEYLEVERETT, LINEARADVECTION, KPP, KPPX, KPPY}; enum TimeIntegrationName{SSPRK2, SSPRK3, SSPRK5, FE}; enum limiterName{NOLIMITER, EXTREMAPRES}; enum ReconstrName{NORECONSTR, ENO, WENO}; enum InterpBaseName{MULTIQUADS, POLYNOMIALS, GAUSSIANS, POLYHARMONICS, NOTUSED}; enum PolyBasis{MONOMIALS, LAGRANGE}; enum FluxName{LAXFRIEDRICHS}; enum QuadratureName{GAUSSLOBATTO}; enum initalCondName{SIN,EULERSOD, BUBBLE, InitKPP, InitKPPY, COSCOS, InitBURGERS, SHOCKX, SINSIN, InitISENTROPICVORTEX - , InitEULERBLAST, InitVORTEXSHOCK, RP, InitDOUBLEMACHREF, InitCUBE, InitFAC, InitAIRFOIL, InitNOZZLE, InitAXISYMMNOZZLE}; -enum sourceName{NOSOURCE, AXISYMMEULER}; + , InitEULERBLAST, InitVORTEXSHOCK, RP, InitDOUBLEMACHREF, InitCUBE, InitFAC, InitAIRFOIL, InitNOZZLE, InitAXISYMMNOZZLE, InitAXISYMMNOZZLERU}; +enum sourceName{NOSOURCE, AXISYMMEULER, AXISYMMEULERRU}; enum stencilChoice{STD, OTF, SYMM, SYMMb, SYMMc}; enum MeshType{TRI, QUAD, HYBRID, TRIQUAD, MIXED, QUADQUAD, RT, RQ}; class Configuration { ModelName model; MeshType meshType; evaluationMethod evalMethod; bcName bcType[6]; ReconstrName reconstr; InterpBaseName interpBase; TimeIntegrationName time_scheme; FluxName fluxScheme; limiterName limiter_name; bcName bcname; initalCondName initialCond; sourceName source; stencilChoice stencilchoice; bool reconstrOnCharact = 0; int order; int degPoly; double tMax; double cfl; double shapeParam; int maxNumWrite = 20; 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; bool load_grid; bool import_sol; std::string grid_file; std::string mesh_type; std::string ic_file; std::string source_file; std::string OUTPUT_PATH; std::string IMPORT_PATH; std::string SOLN_DIR; bool ic_file_read; bool exact_available; int nThreads; std::map constants; public: Configuration(); Configuration(char* file); Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, double cfl); Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, int order, double tMax, std::vector tInclude, double cfl); Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, double cfl); Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, int order, double tMax, std::vector tInclude, double cfl); Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, double cfl); Configuration(ModelName model, evaluationMethod evalMethod, bcName bc, ReconstrName reconstr, InterpBaseName interpBase, bool reconstrOnCharact, int order, double tMax, std::vector tInclude, double cfl); ~Configuration(){}; ModelName getModel(); MeshType getMeshType(); TimeIntegrationName getTimeIntegration(); int getOrder(); evaluationMethod getEvalMethod(); bcName getBc(); ReconstrName getReconstr(); InterpBaseName getInterpBase(); limiterName getLimiter(); FluxName getFlux(); sourceName getSource(); initalCondName getInitialCond(); std::vector getTInclude(); double getTmax(); double getCfl(); int getMaxNumWrite(); bool getReconstrOnCharact(); int getDim(); double getXmax(); double getXmin(); double getYmax(); double getYmin(); double getShapeParam(); std::vector getNx(); std::vector getNy(); void read (); std::string getGridFile(); std::string getOutputPath(); std::string getImportPath(); bool getImportSol(); bool getLoadGrid(); int getNthreads(); stencilChoice getStencilChoice(); std::map getConstants(); int getDegPoly(); 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 read_import (Reader &); }; #endif /* Configuration_hpp */ diff --git a/ConSol/ConSol/Base/Setup.cpp b/ConSol/ConSol/Base/Setup.cpp index c9c7493..175d3e4 100644 --- a/ConSol/ConSol/Base/Setup.cpp +++ b/ConSol/ConSol/Base/Setup.cpp @@ -1,471 +1,478 @@ // // 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 if (config.getMeshType() == HYBRID){ reconstr = new class HybridRBFWENO(); }else if (config.getStencilChoice() == SYMM) { reconstr = new class MRWENO(); }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(); } 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; } 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 { 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; reconstr->initialize(model, interpBase, limiter, mesh, 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; if (source->getIsUsed()){ source->initialize(config); std::cout << "Source term 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 1f92a15..fc52faf 100644 --- a/ConSol/ConSol/Base/Setup.hpp +++ b/ConSol/ConSol/Base/Setup.hpp @@ -1,212 +1,214 @@ // // 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.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 "../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 05f2aaf..032f245 100644 --- a/ConSol/ConSol/CMakeLists.txt +++ b/ConSol/ConSol/CMakeLists.txt @@ -1,316 +1,316 @@ 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 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) + 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) 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/Source/AxisymmEuler_rU_noSwirl.cpp b/ConSol/ConSol/Source/AxisymmEuler_rU_noSwirl.cpp new file mode 100644 index 0000000..c97906a --- /dev/null +++ b/ConSol/ConSol/Source/AxisymmEuler_rU_noSwirl.cpp @@ -0,0 +1,20 @@ +// +// Created by Fabian Moenkeberg on 2020-02-14. +// + +#include "AxisymmEuler_rU_noSwirl.hpp" +std::vector AxisymmEuler_rU_noSwirl::source(std::vector> &U, MeshBase *mesh, double t, double dt, int iCell){ + std::vector ret(U.size(),0); + + long double rrho, rp, rm1, rm2, rE; + rrho = U[iCell][0]; + rm1 = U[iCell][1]; + rm2 = U[iCell][2]; + rE = U[iCell][3]; + rp = (gamma - 1)*(rE - 0.5*(rm1*rm1 + rm2*rm2)/rrho); + + std::vector r = mesh->MidPoint(iCell); + ret[2] = rp/r[1]; + + return ret; +} \ No newline at end of file diff --git a/ConSol/ConSol/Source/AxisymmEuler_rU_noSwirl.hpp b/ConSol/ConSol/Source/AxisymmEuler_rU_noSwirl.hpp new file mode 100644 index 0000000..1e2dd3f --- /dev/null +++ b/ConSol/ConSol/Source/AxisymmEuler_rU_noSwirl.hpp @@ -0,0 +1,33 @@ +// +// Created by Fabian Moenkeberg on 2020-02-14. +// + +#ifndef CONSOL_AXISYMMEULER_RU_NOSWIRL_HPP +#define CONSOL_AXISYMMEULER_RU_NOSWIRL_HPP + + +#include "Source.hpp" + +class AxisymmEuler_rU_noSwirl:public Source { +private: + double gamma = 0; + +public: + AxisymmEuler_rU_noSwirl(){ + this->isUsed = 1; + this->gamma = 1.4; + }; + + AxisymmEuler_rU_noSwirl(ModelBase model){ + this->model = model; + this->isUsed = 1; + this->gamma = 1.4; + }; + + ~AxisymmEuler_rU_noSwirl(){}; + + std::vector source(std::vector> &U, MeshBase *mesh, double t, double dt, int iCell); +}; + + +#endif //CONSOL_AXISYMMEULER_RU_NOSWIRL_HPP diff --git a/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp b/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp index edfe194..b8c9f76 100644 --- a/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp +++ b/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp @@ -1,64 +1,62 @@ // // Created by Fabian Moenkeberg on 2020-02-13. // #include "InitAxiSymmNozzle.hpp" std::vector initAxiSymmNozzle::F(std::vector x){ std::vector ret(4,0); // double r2 = pow(x[0]-xc[0],2) + pow(x[1]-xc[1],2); // ret[0] = pow(1-beta*beta*(gamma-1)*exp(1-r2)/(8*gamma*M_PI*M_PI),1/(gamma-1)); // ret[1] = (M*cos(alpha)-beta*(x[1]-xc[1])*exp((1-r2)*0.5)*0.5/M_PI)*ret[0]; // ret[2] = (M*sin(alpha)+beta*(x[0]-xc[0])*exp((1-r2)*0.5)*0.5/M_PI)*ret[0]; // ret[3] = pow(ret[0],gamma)/(gamma-1)+0.5*(ret[1]*ret[1]+ret[2]*ret[2])/ret[0]; // double r_in = 1; // double u1_in = std::sqrt(gamma)*M; // double u2_in =0; // double p_in = 1; std::vector U_in(4); std::vector U_out(4); U_in.resize(4); U_out.resize(4); double T_in = 300; double NPR = 2.1*3; double m_in = 0.7; double p_in = 101325*NPR; double r_in = p_in/287.14/T_in; double u1_in = m_in/r_in/0.00116; u1_in = 0; double u2_in =0; double E_in = p_in/(gamma - 1) + 0.5*(u1_in*u1_in+u2_in*u2_in)*r_in; double T_out = 300; double u1_out = 0; double u2_out =0; double p_out = 101325; double r_out = 1.3; r_out = p_out/287.14/T_out; double E_out = p_out/(gamma - 1) + 0.5*(u1_out*u1_out+u2_out*u2_out)*r_out; U_in = {r_in, r_in*u1_in, r_in*u2_in, E_in}; U_out = {r_out, r_out*u1_out, r_out*u2_out, E_out}; -// if (false){ -// if (true){ if (x[0] < 110 && x[1]< 40){ - ret[0] = r_in*x[1]; - ret[1] = u1_in*r_in*x[1]; - ret[2] = u2_in*r_in*x[1]; - ret[3] = E_in*x[1]; + ret[0] = r_in; + ret[1] = u1_in*r_in; + ret[2] = u2_in*r_in; + ret[3] = E_in; }else{ - ret[0] = r_out*x[1]; - ret[1] = u1_out*r_out*x[1]; - ret[2] = u2_out*r_out*x[1]; - ret[3] = E_out*x[1]; + ret[0] = r_out; + ret[1] = u1_out*r_out; + ret[2] = u2_out*r_out; + ret[3] = E_out; } return ret; } diff --git a/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp b/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.cpp similarity index 91% copy from ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp copy to ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.cpp index edfe194..c234f51 100644 --- a/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle.cpp +++ b/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.cpp @@ -1,64 +1,65 @@ // -// Created by Fabian Moenkeberg on 2020-02-13. +// Created by Fabian Moenkeberg on 2020-02-14. // -#include "InitAxiSymmNozzle.hpp" -std::vector initAxiSymmNozzle::F(std::vector x){ +#include "InitAxiSymmNozzle_rU.hpp" + +std::vector initAxiSymmNozzle_rU::F(std::vector x){ std::vector ret(4,0); // double r2 = pow(x[0]-xc[0],2) + pow(x[1]-xc[1],2); // ret[0] = pow(1-beta*beta*(gamma-1)*exp(1-r2)/(8*gamma*M_PI*M_PI),1/(gamma-1)); // ret[1] = (M*cos(alpha)-beta*(x[1]-xc[1])*exp((1-r2)*0.5)*0.5/M_PI)*ret[0]; // ret[2] = (M*sin(alpha)+beta*(x[0]-xc[0])*exp((1-r2)*0.5)*0.5/M_PI)*ret[0]; // ret[3] = pow(ret[0],gamma)/(gamma-1)+0.5*(ret[1]*ret[1]+ret[2]*ret[2])/ret[0]; // double r_in = 1; // double u1_in = std::sqrt(gamma)*M; // double u2_in =0; // double p_in = 1; std::vector U_in(4); std::vector U_out(4); U_in.resize(4); U_out.resize(4); double T_in = 300; double NPR = 2.1*3; double m_in = 0.7; double p_in = 101325*NPR; double r_in = p_in/287.14/T_in; double u1_in = m_in/r_in/0.00116; u1_in = 0; double u2_in =0; double E_in = p_in/(gamma - 1) + 0.5*(u1_in*u1_in+u2_in*u2_in)*r_in; double T_out = 300; double u1_out = 0; double u2_out =0; double p_out = 101325; double r_out = 1.3; r_out = p_out/287.14/T_out; double E_out = p_out/(gamma - 1) + 0.5*(u1_out*u1_out+u2_out*u2_out)*r_out; U_in = {r_in, r_in*u1_in, r_in*u2_in, E_in}; U_out = {r_out, r_out*u1_out, r_out*u2_out, E_out}; // if (false){ // if (true){ if (x[0] < 110 && x[1]< 40){ ret[0] = r_in*x[1]; ret[1] = u1_in*r_in*x[1]; ret[2] = u2_in*r_in*x[1]; ret[3] = E_in*x[1]; }else{ ret[0] = r_out*x[1]; ret[1] = u1_out*r_out*x[1]; ret[2] = u2_out*r_out*x[1]; ret[3] = E_out*x[1]; } return ret; -} +} \ No newline at end of file diff --git a/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.hpp b/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.hpp new file mode 100644 index 0000000..97b9931 --- /dev/null +++ b/ConSol/ConSol/Utility/InitialConditions/2D/InitAxiSymmNozzle_rU.hpp @@ -0,0 +1,46 @@ +// +// Created by Fabian Moenkeberg on 2020-02-14. +// + +#ifndef CONSOL_INITAXISYMMNOZZLE_RU_HPP +#define CONSOL_INITAXISYMMNOZZLE_RU_HPP + +#include +#include "../../AbstractFunct.hpp" + +class initAxiSymmNozzle_rU :public AbstractFunct { + std::vector xc; + double alpha; + double beta; + double gamma; + double M; +public: + initAxiSymmNozzle_rU(){ + dim = 2; + xc.resize(2); + xc[0] = 0; + xc[1] = 0; + alpha = 0; + M = 0.5; + gamma = 1.4; + }; + + initAxiSymmNozzle_rU(double alpha, double M){ + dim = 2; + xc.resize(2); + xc[0] = 0; + xc[1] = 0; + this->alpha = alpha; + this->beta = beta; + this->M = M; + this->gamma = 1.4; + }; + + ~initAxiSymmNozzle_rU(){}; + + std::vector F(std::vector x); +}; + + + +#endif //CONSOL_INITAXISYMMNOZZLE_RU_HPP