1 /// \file Sedov.cpp 2 /// \brief Sedov blast test 3 /// \author Pavel Sevecek (sevecek at sirrah.troja.mff.cuni.cz) 4 /// \date 2016-2021 5 6 #include "Sph.h" 7 #include "catch.hpp" 8 #include <iostream> 9 10 using namespace Sph; 11 12 class PlanarDistribution : public IDistribution { 13 public: generate(IScheduler &,const Size n,const IDomain & domain) const14 virtual Array<Vector> generate(IScheduler&, const Size n, const IDomain& domain) const override { 15 const Box box = domain.getBoundingBox(); 16 const Float h = sqrt(1._f / n); 17 const Float dx = h; 18 const Float dy = sqrt(3._f) / 2._f * dx; 19 Array<Vector> r; 20 int row = 0; 21 for (Float y = box.lower()[Y]; y <= box.upper()[Y]; y += dy, row++) { 22 for (Float x = box.lower()[X]; x <= box.upper()[X]; x += dx) { 23 Float delta = (row % 2) ? 0.5_f * dx : 0._f; 24 r.push(Vector(x + delta, y, 0._f, h)); 25 } 26 } 27 return r; 28 } 29 }; 30 31 class Sedov : public IRun { 32 public: Sedov()33 Sedov() { 34 // Global settings of the problem 35 this->settings.set(RunSettingsId::RUN_NAME, String("Sedov Blast Problem")) 36 .set(RunSettingsId::RUN_END_TIME, 8._f) 37 .set(RunSettingsId::RUN_OUTPUT_TYPE, IoEnum::TEXT_FILE) 38 .set(RunSettingsId::RUN_OUTPUT_INTERVAL, 0.08_f) 39 .set(RunSettingsId::RUN_OUTPUT_PATH, String("")) 40 .set(RunSettingsId::RUN_OUTPUT_NAME, String("sedov/sedov_%d.txt")) 41 .set(RunSettingsId::SPH_AV_ALPHA, 1.5_f) 42 .set(RunSettingsId::SPH_AV_BETA, 3._f) 43 .set(RunSettingsId::SPH_KERNEL, KernelEnum::CUBIC_SPLINE) 44 .set(RunSettingsId::TIMESTEPPING_INTEGRATOR, TimesteppingEnum::EULER_EXPLICIT) 45 .set(RunSettingsId::TIMESTEPPING_INITIAL_TIMESTEP, 1.e-5_f) 46 .set(RunSettingsId::TIMESTEPPING_MAX_TIMESTEP, 0.1_f) 47 .set(RunSettingsId::TIMESTEPPING_COURANT_NUMBER, 0.2_f) 48 .set(RunSettingsId::TIMESTEPPING_CRITERION, TimeStepCriterionEnum::COURANT) 49 .set(RunSettingsId::SPH_SOLVER_FORCES, ForceEnum::PRESSURE) 50 .set(RunSettingsId::SPH_USE_AC, true) 51 .set(RunSettingsId::SPH_FINDER, FinderEnum::UNIFORM_GRID) 52 .set(RunSettingsId::SPH_ADAPTIVE_SMOOTHING_LENGTH, SmoothingLengthEnum::CONTINUITY_EQUATION) 53 .set(RunSettingsId::DOMAIN_TYPE, DomainEnum::BLOCK) 54 //.set(RunSettingsId::DOMAIN_BOUNDARY, BoundaryEnum::FROZEN_PARTICLES) 55 .set(RunSettingsId::DOMAIN_SIZE, Vector(1._f)); 56 } 57 setUp(SharedPtr<Storage> storage)58 virtual void setUp(SharedPtr<Storage> storage) override { 59 BodySettings body; 60 body.set(BodySettingsId::DENSITY, 1._f) 61 .set(BodySettingsId::DENSITY_RANGE, Interval(1.e-3_f, INFTY)) 62 .set(BodySettingsId::ENERGY, 0._f) 63 .set(BodySettingsId::ENERGY_RANGE, Interval(0._f, INFTY)) 64 .set(BodySettingsId::EOS, EosEnum::IDEAL_GAS) 65 .set(BodySettingsId::RHEOLOGY_YIELDING, YieldingEnum::NONE) 66 .set(BodySettingsId::RHEOLOGY_DAMAGE, FractureEnum::NONE) 67 .set(BodySettingsId::ADIABATIC_INDEX, 1.6666666666_f) 68 .set(BodySettingsId::DISTRIBUTE_MODE_SPH5, true); 69 70 *storage = Storage(Factory::getMaterial(body)); 71 AutoPtr<IDistribution> dist = makeAuto<PlanarDistribution>(); // Factory::getDistribution(body); 72 BlockDomain domain(Vector(0._f), Vector(1._f, 1._f, 1.e-3_f)); 73 Array<Vector> pos = dist->generate(*scheduler, 100000, domain); 74 const Float eta = 1.5_f; 75 for (Size i = 0; i < pos.size(); ++i) { 76 pos[i][H] *= eta; 77 } 78 const Float m = body.get<Float>(BodySettingsId::DENSITY) / pos.size(); 79 80 storage->insert<Vector>(QuantityId::POSITION, OrderEnum::SECOND, std::move(pos)); 81 storage->insert<Float>(QuantityId::MASS, OrderEnum::ZERO, 82 m); // rho * S / N 83 84 EquationHolder eqs = getStandardEquations(settings); 85 this->solver = makeAuto<SymmetricSolver<2>>(*scheduler, settings, eqs); 86 IMaterial& mat = storage->getMaterial(0); 87 solver->create(*storage, mat); 88 MaterialInitialContext context(settings); 89 mat.create(*storage, context); 90 91 ArrayView<Float> u = storage->getValue<Float>(QuantityId::ENERGY); 92 ArrayView<Vector> r = storage->getValue<Vector>(QuantityId::POSITION); 93 Float E0 = 0._f; 94 for (Size i = 0; i < r.size(); ++i) { 95 if (getLength(r[i]) < 0.015_f) { 96 u[i] = 4._f; 97 E0 += u[i] * m; 98 } 99 } 100 std::cout << "E0 = " << E0 << std::endl; 101 } 102 103 protected: tearDown(const Storage & UNUSED (storage),const Statistics & UNUSED (stats))104 virtual void tearDown(const Storage& UNUSED(storage), const Statistics& UNUSED(stats)) override {} 105 }; 106 107 TEST_CASE("Sedov", "[sedov]") { 108 Sedov run; 109 Storage storage; 110 run.run(storage); 111 } 112