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