1 #include "run/jobs/SimulationJobs.h"
2 #include "gravity/AggregateSolver.h"
3 #include "io/LogWriter.h"
4 #include "io/Logger.h"
5 #include "io/Output.h"
6 #include "run/IRun.h"
7 #include "run/SpecialEntries.h"
8 #include "sph/solvers/StabilizationSolver.h"
9
10 NAMESPACE_SPH_BEGIN
11
12 /// \todo generailize, add generic triggers to UI
13 class EnergyLogWriter : public ILogWriter {
14 public:
15 using ILogWriter::ILogWriter;
16
17 private:
write(const Storage & storage,const Statistics & stats)18 virtual void write(const Storage& storage, const Statistics& stats) override {
19 const Float t = stats.get<Float>(StatisticsId::RUN_TIME);
20 const Float e = TotalEnergy().evaluate(storage);
21 logger->write(t, " ", e);
22 }
23 };
24
getIdentifier(const String & name)25 static String getIdentifier(const String& name) {
26 String escaped = name;
27 escaped.replaceAll(" ", "-");
28 return escaped.toLowercase();
29 }
30
31 // ----------------------------------------------------------------------------------------------------------
32 // SphJob
33 // ----------------------------------------------------------------------------------------------------------
34
overrideSettings(const RunSettings & settings,const RunSettings & overrides,const bool isResumed)35 static RunSettings overrideSettings(const RunSettings& settings,
36 const RunSettings& overrides,
37 const bool isResumed) {
38 RunSettings actual = settings;
39 actual.addEntries(overrides);
40
41 if (!isResumed) {
42 // reset the (potentially) overriden values back to original
43 actual.set(RunSettingsId::RUN_START_TIME, settings.get<Float>(RunSettingsId::RUN_START_TIME));
44 actual.set(RunSettingsId::TIMESTEPPING_INITIAL_TIMESTEP,
45 settings.get<Float>(RunSettingsId::TIMESTEPPING_INITIAL_TIMESTEP));
46 actual.set(
47 RunSettingsId::RUN_OUTPUT_FIRST_INDEX, settings.get<int>(RunSettingsId::RUN_OUTPUT_FIRST_INDEX));
48 }
49 return actual;
50 }
51
addTimeSteppingCategory(VirtualSettings & connector,RunSettings & settings,bool & resumeRun)52 static void addTimeSteppingCategory(VirtualSettings& connector, RunSettings& settings, bool& resumeRun) {
53 auto courantEnabler = [&settings] {
54 Flags<TimeStepCriterionEnum> criteria =
55 settings.getFlags<TimeStepCriterionEnum>(RunSettingsId::TIMESTEPPING_CRITERION);
56 return criteria.has(TimeStepCriterionEnum::COURANT);
57 };
58 auto derivativeEnabler = [&settings] {
59 Flags<TimeStepCriterionEnum> criteria =
60 settings.getFlags<TimeStepCriterionEnum>(RunSettingsId::TIMESTEPPING_CRITERION);
61 return criteria.hasAny(TimeStepCriterionEnum::DERIVATIVES, TimeStepCriterionEnum::ACCELERATION);
62 };
63 auto divergenceEnabler = [&settings] {
64 Flags<TimeStepCriterionEnum> criteria =
65 settings.getFlags<TimeStepCriterionEnum>(RunSettingsId::TIMESTEPPING_CRITERION);
66 return criteria.has(TimeStepCriterionEnum::DIVERGENCE);
67 };
68
69 VirtualSettings::Category& rangeCat = connector.addCategory("Integration");
70 rangeCat.connect<Float>("Duration [s]", settings, RunSettingsId::RUN_END_TIME);
71 rangeCat.connect("Use start time of input", "is_resumed", resumeRun);
72 rangeCat.connect<Float>("Maximal timestep [s]", settings, RunSettingsId::TIMESTEPPING_MAX_TIMESTEP);
73 rangeCat.connect<Float>("Initial timestep [s]", settings, RunSettingsId::TIMESTEPPING_INITIAL_TIMESTEP);
74 rangeCat.connect<EnumWrapper>("Integrator", settings, RunSettingsId::TIMESTEPPING_INTEGRATOR);
75 rangeCat.connect<Flags<TimeStepCriterionEnum>>(
76 "Time step criteria", settings, RunSettingsId::TIMESTEPPING_CRITERION);
77 rangeCat.connect<Float>("Courant number", settings, RunSettingsId::TIMESTEPPING_COURANT_NUMBER)
78 .setEnabler(courantEnabler);
79 rangeCat.connect<Float>("Derivative factor", settings, RunSettingsId::TIMESTEPPING_DERIVATIVE_FACTOR)
80 .setEnabler(derivativeEnabler);
81 rangeCat.connect<Float>("Divergence factor", settings, RunSettingsId::TIMESTEPPING_DIVERGENCE_FACTOR)
82 .setEnabler(divergenceEnabler);
83 rangeCat.connect<Float>("Max step change", settings, RunSettingsId::TIMESTEPPING_MAX_INCREASE);
84 }
85
addGravityCategory(VirtualSettings & connector,RunSettings & settings)86 static void addGravityCategory(VirtualSettings& connector, RunSettings& settings) {
87 VirtualSettings::Category& gravityCat = connector.addCategory("Gravity");
88 gravityCat.connect<EnumWrapper>("Gravity solver", settings, RunSettingsId::GRAVITY_SOLVER);
89 gravityCat.connect<Float>("Opening angle", settings, RunSettingsId::GRAVITY_OPENING_ANGLE)
90 .setEnabler([&settings] {
91 return settings.get<GravityEnum>(RunSettingsId::GRAVITY_SOLVER) == GravityEnum::BARNES_HUT;
92 });
93 gravityCat.connect<int>("Multipole order", settings, RunSettingsId::GRAVITY_MULTIPOLE_ORDER);
94 gravityCat.connect<EnumWrapper>("Softening kernel", settings, RunSettingsId::GRAVITY_KERNEL);
95 gravityCat.connect<Float>(
96 "Recomputation period [s]", settings, RunSettingsId::GRAVITY_RECOMPUTATION_PERIOD);
97 }
98
addOutputCategory(VirtualSettings & connector,RunSettings & settings,const SharedToken & owner)99 static void addOutputCategory(VirtualSettings& connector, RunSettings& settings, const SharedToken& owner) {
100 auto enabler = [&settings] {
101 const IoEnum type = settings.get<IoEnum>(RunSettingsId::RUN_OUTPUT_TYPE);
102 return type != IoEnum::NONE;
103 };
104
105 VirtualSettings::Category& outputCat = connector.addCategory("Output");
106 outputCat.connect<EnumWrapper>("Format", settings, RunSettingsId::RUN_OUTPUT_TYPE)
107 .setValidator([](const IVirtualEntry::Value& value) {
108 const IoEnum type = IoEnum(value.get<EnumWrapper>());
109 return type == IoEnum::NONE || getIoCapabilities(type).has(IoCapability::OUTPUT);
110 })
111 .addAccessor(owner,
112 [&settings](const IVirtualEntry::Value& value) {
113 const IoEnum type = IoEnum(value.get<EnumWrapper>());
114 Path name = Path(settings.get<String>(RunSettingsId::RUN_OUTPUT_NAME));
115 if (Optional<String> extension = getIoExtension(type)) {
116 name.replaceExtension(extension.value());
117 }
118 settings.set(RunSettingsId::RUN_OUTPUT_NAME, name.string());
119 })
120 .setSideEffect(); // needs to update the 'File mask' entry
121 outputCat.connect<Path>("Directory", settings, RunSettingsId::RUN_OUTPUT_PATH)
122 .setEnabler(enabler)
123 .setPathType(IVirtualEntry::PathType::DIRECTORY);
124 outputCat.connect<String>("File mask", settings, RunSettingsId::RUN_OUTPUT_NAME).setEnabler(enabler);
125 outputCat.connect<Flags<OutputQuantityFlag>>("Quantities", settings, RunSettingsId::RUN_OUTPUT_QUANTITIES)
126 .setEnabler([&settings] {
127 const IoEnum type = settings.get<IoEnum>(RunSettingsId::RUN_OUTPUT_TYPE);
128 return type == IoEnum::TEXT_FILE || type == IoEnum::VTK_FILE;
129 });
130 outputCat.connect<EnumWrapper>("Output spacing", settings, RunSettingsId::RUN_OUTPUT_SPACING)
131 .setEnabler(enabler);
132 outputCat.connect<Float>("Output interval [s]", settings, RunSettingsId::RUN_OUTPUT_INTERVAL)
133 .setEnabler([&] {
134 const IoEnum type = settings.get<IoEnum>(RunSettingsId::RUN_OUTPUT_TYPE);
135 const OutputSpacing spacing = settings.get<OutputSpacing>(RunSettingsId::RUN_OUTPUT_SPACING);
136 return type != IoEnum::NONE && spacing != OutputSpacing::CUSTOM;
137 });
138 outputCat.connect<String>("Custom times [s]", settings, RunSettingsId::RUN_OUTPUT_CUSTOM_TIMES)
139 .setEnabler([&] {
140 const IoEnum type = settings.get<IoEnum>(RunSettingsId::RUN_OUTPUT_TYPE);
141 const OutputSpacing spacing = settings.get<OutputSpacing>(RunSettingsId::RUN_OUTPUT_SPACING);
142 return type != IoEnum::NONE && spacing == OutputSpacing::CUSTOM;
143 });
144 }
145
addLoggerCategory(VirtualSettings & connector,RunSettings & settings)146 static void addLoggerCategory(VirtualSettings& connector, RunSettings& settings) {
147 VirtualSettings::Category& loggerCat = connector.addCategory("Logging");
148 loggerCat.connect<EnumWrapper>("Logger", settings, RunSettingsId::RUN_LOGGER);
149 loggerCat.connect<Path>("Log file", settings, RunSettingsId::RUN_LOGGER_FILE)
150 .setPathType(IVirtualEntry::PathType::OUTPUT_FILE)
151 .setEnabler(
152 [&settings] { return settings.get<LoggerEnum>(RunSettingsId::RUN_LOGGER) == LoggerEnum::FILE; });
153 loggerCat.connect<int>("Log verbosity", settings, RunSettingsId::RUN_LOGGER_VERBOSITY);
154 }
155
156
157 class SphRun : public IRun {
158 protected:
159 SharedPtr<IDomain> domain;
160
161 public:
SphRun(const RunSettings & run,SharedPtr<IDomain> domain)162 explicit SphRun(const RunSettings& run, SharedPtr<IDomain> domain)
163 : domain(domain) {
164 settings = run;
165 scheduler = Factory::getScheduler(settings);
166 }
167
setUp(SharedPtr<Storage> storage)168 virtual void setUp(SharedPtr<Storage> storage) override {
169 AutoPtr<IBoundaryCondition> bc = Factory::getBoundaryConditions(settings, domain);
170 solver = Factory::getSolver(*scheduler, settings, std::move(bc));
171
172 for (Size matId = 0; matId < storage->getMaterialCnt(); ++matId) {
173 solver->create(*storage, storage->getMaterial(matId));
174 }
175 }
176
tearDown(const Storage & storage,const Statistics & stats)177 virtual void tearDown(const Storage& storage, const Statistics& stats) override {
178 // last dump after simulation ends
179 output->dump(storage, stats);
180 }
181 };
182
SphJob(const String & name,const RunSettings & overrides)183 SphJob::SphJob(const String& name, const RunSettings& overrides)
184 : IRunJob(name) {
185 settings = getDefaultSettings(name);
186
187 settings.addEntries(overrides);
188 }
189
getDefaultSettings(const String & name)190 RunSettings SphJob::getDefaultSettings(const String& name) {
191 const Size dumpCnt = 10;
192 const Interval timeRange(0, 10);
193
194 RunSettings settings;
195 settings.set(RunSettingsId::TIMESTEPPING_INTEGRATOR, TimesteppingEnum::PREDICTOR_CORRECTOR)
196 .set(RunSettingsId::TIMESTEPPING_INITIAL_TIMESTEP, 0.01_f)
197 .set(RunSettingsId::TIMESTEPPING_MAX_TIMESTEP, 10._f)
198 .set(RunSettingsId::TIMESTEPPING_COURANT_NUMBER, 0.2_f)
199 .set(RunSettingsId::RUN_START_TIME, timeRange.lower())
200 .set(RunSettingsId::RUN_END_TIME, timeRange.upper())
201 .set(RunSettingsId::RUN_NAME, name)
202 .set(RunSettingsId::RUN_OUTPUT_INTERVAL, timeRange.size() / dumpCnt)
203 .set(RunSettingsId::RUN_OUTPUT_TYPE, IoEnum::NONE)
204 .set(RunSettingsId::RUN_OUTPUT_NAME, getIdentifier(name) + "_%d.ssf")
205 .set(RunSettingsId::RUN_VERBOSE_NAME, getIdentifier(name) + ".log")
206 .set(RunSettingsId::SPH_SOLVER_TYPE, SolverEnum::ASYMMETRIC_SOLVER)
207 .set(RunSettingsId::SPH_SOLVER_FORCES,
208 ForceEnum::PRESSURE | ForceEnum::SOLID_STRESS | ForceEnum::SELF_GRAVITY)
209 .set(RunSettingsId::SPH_DISCRETIZATION, DiscretizationEnum::STANDARD)
210 .set(RunSettingsId::SPH_FINDER, FinderEnum::KD_TREE)
211 .set(RunSettingsId::SPH_AV_TYPE, ArtificialViscosityEnum::STANDARD)
212 .set(RunSettingsId::SPH_AV_ALPHA, 1.5_f)
213 .set(RunSettingsId::SPH_AV_BETA, 3._f)
214 .set(RunSettingsId::SPH_KERNEL, KernelEnum::CUBIC_SPLINE)
215 .set(RunSettingsId::GRAVITY_SOLVER, GravityEnum::BARNES_HUT)
216 .set(RunSettingsId::GRAVITY_KERNEL, GravityKernelEnum::SPH_KERNEL)
217 .set(RunSettingsId::GRAVITY_OPENING_ANGLE, 0.8_f)
218 .set(RunSettingsId::GRAVITY_RECOMPUTATION_PERIOD, 5._f)
219 .set(RunSettingsId::FINDER_LEAF_SIZE, 20)
220 .set(RunSettingsId::SPH_STABILIZATION_DAMPING, 0.1_f)
221 .set(RunSettingsId::RUN_THREAD_GRANULARITY, 1000)
222 .set(RunSettingsId::SPH_ADAPTIVE_SMOOTHING_LENGTH, EMPTY_FLAGS)
223 .set(RunSettingsId::SPH_ASYMMETRIC_COMPUTE_RADII_HASH_MAP, false)
224 .set(RunSettingsId::SPH_STRAIN_RATE_CORRECTION_TENSOR, true)
225 .set(RunSettingsId::RUN_DIAGNOSTICS_INTERVAL, 1._f);
226 return settings;
227 }
228
getSettings()229 VirtualSettings SphJob::getSettings() {
230 VirtualSettings connector;
231 addGenericCategory(connector, instName);
232 addTimeSteppingCategory(connector, settings, isResumed);
233
234 auto stressEnabler = [this] {
235 return settings.getFlags<ForceEnum>(RunSettingsId::SPH_SOLVER_FORCES).has(ForceEnum::SOLID_STRESS);
236 };
237 auto avEnabler = [this] {
238 return settings.get<ArtificialViscosityEnum>(RunSettingsId::SPH_AV_TYPE) !=
239 ArtificialViscosityEnum::NONE;
240 };
241 auto asEnabler = [this] { return settings.get<bool>(RunSettingsId::SPH_AV_USE_STRESS); };
242 // auto acEnabler = [this] { return settings.get<bool>(RunSettingsId::SPH_USE_AC); };
243 auto deltaSphEnabler = [this] { return settings.get<bool>(RunSettingsId::SPH_USE_DELTASPH); };
244 auto enforceEnabler = [this] {
245 return settings.getFlags<SmoothingLengthEnum>(RunSettingsId::SPH_ADAPTIVE_SMOOTHING_LENGTH)
246 .has(SmoothingLengthEnum::SOUND_SPEED_ENFORCING);
247 };
248
249 VirtualSettings::Category& solverCat = connector.addCategory("SPH solver");
250 solverCat.connect<Flags<ForceEnum>>("Forces", settings, RunSettingsId::SPH_SOLVER_FORCES);
251 solverCat.connect<Vector>("Constant acceleration", settings, RunSettingsId::FRAME_CONSTANT_ACCELERATION);
252 solverCat.connect<Float>("Tides mass [M_earth]", settings, RunSettingsId::FRAME_TIDES_MASS)
253 .setUnits(Constants::M_earth);
254 solverCat.connect<Vector>("Tides position [R_earth]", settings, RunSettingsId::FRAME_TIDES_POSITION)
255 .setUnits(Constants::R_earth);
256 solverCat.connect<EnumWrapper>("Solver type", settings, RunSettingsId::SPH_SOLVER_TYPE);
257 solverCat.connect<EnumWrapper>("SPH discretization", settings, RunSettingsId::SPH_DISCRETIZATION);
258 solverCat.connect<Flags<SmoothingLengthEnum>>(
259 "Adaptive smoothing length", settings, RunSettingsId::SPH_ADAPTIVE_SMOOTHING_LENGTH);
260 solverCat.connect<Float>("Minimal smoothing length", settings, RunSettingsId::SPH_SMOOTHING_LENGTH_MIN)
261 .setEnabler([this] {
262 return settings.getFlags<SmoothingLengthEnum>(RunSettingsId::SPH_ADAPTIVE_SMOOTHING_LENGTH) !=
263 EMPTY_FLAGS;
264 });
265 solverCat
266 .connect<Float>("Neighbor count enforcing strength", settings, RunSettingsId::SPH_NEIGHBOR_ENFORCING)
267 .setEnabler(enforceEnabler);
268 solverCat.connect<Interval>("Neighbor range", settings, RunSettingsId::SPH_NEIGHBOR_RANGE)
269 .setEnabler(enforceEnabler);
270 solverCat
271 .connect<bool>("Use radii hash map", settings, RunSettingsId::SPH_ASYMMETRIC_COMPUTE_RADII_HASH_MAP)
272 .setEnabler([this] {
273 return settings.get<SolverEnum>(RunSettingsId::SPH_SOLVER_TYPE) == SolverEnum::ASYMMETRIC_SOLVER;
274 });
275 solverCat
276 .connect<bool>("Apply correction tensor", settings, RunSettingsId::SPH_STRAIN_RATE_CORRECTION_TENSOR)
277 .setEnabler(stressEnabler);
278 solverCat.connect<bool>("Sum only undamaged particles", settings, RunSettingsId::SPH_SUM_ONLY_UNDAMAGED);
279 solverCat.connect<EnumWrapper>("Continuity mode", settings, RunSettingsId::SPH_CONTINUITY_MODE);
280 solverCat.connect<EnumWrapper>("Neighbor finder", settings, RunSettingsId::SPH_FINDER);
281 solverCat.connect<EnumWrapper>("Boundary condition", settings, RunSettingsId::DOMAIN_BOUNDARY);
282
283 VirtualSettings::Category& avCat = connector.addCategory("Artificial viscosity");
284 avCat.connect<EnumWrapper>("Artificial viscosity type", settings, RunSettingsId::SPH_AV_TYPE);
285 avCat.connect<bool>("Apply Balsara switch", settings, RunSettingsId::SPH_AV_USE_BALSARA)
286 .setEnabler(avEnabler);
287 avCat.connect<Float>("Artificial viscosity alpha", settings, RunSettingsId::SPH_AV_ALPHA)
288 .setEnabler(avEnabler);
289 avCat.connect<Float>("Artificial viscosity beta", settings, RunSettingsId::SPH_AV_BETA)
290 .setEnabler(avEnabler);
291 avCat.connect<bool>("Apply artificial stress", settings, RunSettingsId::SPH_AV_USE_STRESS);
292 avCat.connect<Float>("Artificial stress factor", settings, RunSettingsId::SPH_AV_STRESS_FACTOR)
293 .setEnabler(asEnabler);
294 avCat.connect<Float>("Artificial stress exponent", settings, RunSettingsId::SPH_AV_STRESS_EXPONENT)
295 .setEnabler(asEnabler);
296 avCat.connect<bool>("Apply artificial conductivity", settings, RunSettingsId::SPH_USE_AC);
297 avCat.connect<EnumWrapper>("Signal speed", settings, RunSettingsId::SPH_AC_SIGNAL_SPEED)
298 .setEnabler([this] { return settings.get<bool>(RunSettingsId::SPH_USE_AC); });
299
300 VirtualSettings::Category& modCat = connector.addCategory("SPH modifications");
301 modCat.connect<bool>("Enable XPSH", settings, RunSettingsId::SPH_USE_XSPH);
302 modCat.connect<Float>("XSPH epsilon", settings, RunSettingsId::SPH_XSPH_EPSILON).setEnabler([this] {
303 return settings.get<bool>(RunSettingsId::SPH_USE_XSPH);
304 });
305 modCat.connect<bool>("Enable delta-SPH", settings, RunSettingsId::SPH_USE_DELTASPH);
306 modCat.connect<Float>("delta-SPH alpha", settings, RunSettingsId::SPH_VELOCITY_DIFFUSION_ALPHA)
307 .setEnabler(deltaSphEnabler);
308 modCat.connect<Float>("delta-SPH delta", settings, RunSettingsId::SPH_DENSITY_DIFFUSION_DELTA)
309 .setEnabler(deltaSphEnabler);
310
311 auto scriptEnabler = [this] { return settings.get<bool>(RunSettingsId::SPH_SCRIPT_ENABLE); };
312
313 VirtualSettings::Category& scriptCat = connector.addCategory("Scripts");
314 scriptCat.connect<bool>("Enable script", settings, RunSettingsId::SPH_SCRIPT_ENABLE);
315 scriptCat.connect<Path>("Script file", settings, RunSettingsId::SPH_SCRIPT_FILE)
316 .setEnabler(scriptEnabler)
317 .setPathType(IVirtualEntry::PathType::INPUT_FILE)
318 .setFileFormats({ { "Chaiscript script", "chai" } });
319 scriptCat.connect<Float>("Script period [s]", settings, RunSettingsId::SPH_SCRIPT_PERIOD)
320 .setEnabler(scriptEnabler);
321 scriptCat.connect<bool>("Run only once", settings, RunSettingsId::SPH_SCRIPT_ONESHOT)
322 .setEnabler(scriptEnabler);
323
324 addGravityCategory(connector, settings);
325 addOutputCategory(connector, settings, *this);
326 addLoggerCategory(connector, settings);
327
328 return connector;
329 }
330
getRun(const RunSettings & overrides) const331 AutoPtr<IRun> SphJob::getRun(const RunSettings& overrides) const {
332 SPH_ASSERT(overrides.size() < 15); // not really required, just checking that we don't override everything
333 const BoundaryEnum boundary = settings.get<BoundaryEnum>(RunSettingsId::DOMAIN_BOUNDARY);
334 SharedPtr<IDomain> domain;
335 if (boundary != BoundaryEnum::NONE) {
336 domain = this->getInput<IDomain>("boundary");
337 }
338
339 RunSettings run = overrideSettings(settings, overrides, isResumed);
340 if (!run.getFlags<ForceEnum>(RunSettingsId::SPH_SOLVER_FORCES).has(ForceEnum::SOLID_STRESS)) {
341 run.set(RunSettingsId::SPH_STRAIN_RATE_CORRECTION_TENSOR, false);
342 }
343
344 return makeAuto<SphRun>(run, domain);
345 }
346
347 static JobRegistrar sRegisterSph(
348 "SPH run",
349 "simulations",
__anon0f2951f41602(const String& name) 350 [](const String& name) { return makeAuto<SphJob>(name, EMPTY_SETTINGS); },
351 "Runs a SPH simulation, using provided initial conditions.");
352
353 // ----------------------------------------------------------------------------------------------------------
354 // SphStabilizationJob
355 // ----------------------------------------------------------------------------------------------------------
356
357 class SphStabilizationRun : public SphRun {
358 public:
359 using SphRun::SphRun;
360
setUp(SharedPtr<Storage> storage)361 virtual void setUp(SharedPtr<Storage> storage) override {
362 AutoPtr<IBoundaryCondition> bc = Factory::getBoundaryConditions(settings, domain);
363 solver = makeAuto<StabilizationSolver>(*scheduler, settings, std::move(bc));
364
365 for (Size matId = 0; matId < storage->getMaterialCnt(); ++matId) {
366 solver->create(*storage, storage->getMaterial(matId));
367 }
368 }
369 };
370
getSettings()371 VirtualSettings SphStabilizationJob::getSettings() {
372 VirtualSettings connector = SphJob::getSettings();
373
374 VirtualSettings::Category& stabCat = connector.addCategory("Stabilization");
375 stabCat.connect<Float>("Damping coefficient", settings, RunSettingsId::SPH_STABILIZATION_DAMPING);
376
377 return connector;
378 }
379
getRun(const RunSettings & overrides) const380 AutoPtr<IRun> SphStabilizationJob::getRun(const RunSettings& overrides) const {
381 RunSettings run = overrideSettings(settings, overrides, isResumed);
382 const BoundaryEnum boundary = settings.get<BoundaryEnum>(RunSettingsId::DOMAIN_BOUNDARY);
383 SharedPtr<IDomain> domain;
384 if (boundary != BoundaryEnum::NONE) {
385 domain = this->getInput<IDomain>("boundary");
386 }
387 return makeAuto<SphStabilizationRun>(run, domain);
388 }
389
390 static JobRegistrar sRegisterSphStab(
391 "SPH stabilization",
392 "stabilization",
393 "simulations",
__anon0f2951f41702(const String& name) 394 [](const String& name) { return makeAuto<SphStabilizationJob>(name, EMPTY_SETTINGS); },
395 "Runs a SPH simulation with a damping term, suitable for stabilization of non-equilibrium initial "
396 "conditions.");
397
398
399 // ----------------------------------------------------------------------------------------------------------
400 // NBodyJob
401 // ----------------------------------------------------------------------------------------------------------
402
403 class NBodyRun : public IRun {
404 private:
405 bool useSoft;
406
407 public:
NBodyRun(const RunSettings & run,const bool useSoft)408 NBodyRun(const RunSettings& run, const bool useSoft)
409 : useSoft(useSoft) {
410 settings = run;
411 scheduler = Factory::getScheduler(settings);
412 }
413
setUp(SharedPtr<Storage> storage)414 virtual void setUp(SharedPtr<Storage> storage) override {
415 logger = Factory::getLogger(settings);
416
417 const bool aggregateEnable = settings.get<bool>(RunSettingsId::NBODY_AGGREGATES_ENABLE);
418 const AggregateEnum aggregateSource =
419 settings.get<AggregateEnum>(RunSettingsId::NBODY_AGGREGATES_SOURCE);
420 if (aggregateEnable) {
421 AutoPtr<AggregateSolver> aggregates = makeAuto<AggregateSolver>(*scheduler, settings);
422 aggregates->createAggregateData(*storage, aggregateSource);
423 solver = std::move(aggregates);
424 } else if (useSoft) {
425 solver = makeAuto<SoftSphereSolver>(*scheduler, settings);
426 } else {
427 solver = makeAuto<HardSphereSolver>(*scheduler, settings);
428 }
429
430 NullMaterial mtl(BodySettings::getDefaults());
431 solver->create(*storage, mtl);
432
433 setPersistentIndices(*storage);
434 }
435
tearDown(const Storage & storage,const Statistics & stats)436 virtual void tearDown(const Storage& storage, const Statistics& stats) override {
437 output->dump(storage, stats);
438 }
439 };
440
NBodyJob(const String & name,const RunSettings & overrides)441 NBodyJob::NBodyJob(const String& name, const RunSettings& overrides)
442 : IRunJob(name) {
443
444 settings = getDefaultSettings(name);
445 settings.addEntries(overrides);
446 }
447
getDefaultSettings(const String & name)448 RunSettings NBodyJob::getDefaultSettings(const String& name) {
449 const Interval timeRange(0, 1.e6_f);
450 RunSettings settings;
451 settings.set(RunSettingsId::RUN_NAME, name)
452 .set(RunSettingsId::RUN_TYPE, RunTypeEnum::NBODY)
453 .set(RunSettingsId::TIMESTEPPING_INTEGRATOR, TimesteppingEnum::LEAP_FROG)
454 .set(RunSettingsId::TIMESTEPPING_INITIAL_TIMESTEP, 0.01_f)
455 .set(RunSettingsId::TIMESTEPPING_MAX_TIMESTEP, 10._f)
456 .set(RunSettingsId::TIMESTEPPING_CRITERION, TimeStepCriterionEnum::ACCELERATION)
457 .set(RunSettingsId::TIMESTEPPING_DERIVATIVE_FACTOR, 0.2_f)
458 .set(RunSettingsId::RUN_START_TIME, timeRange.lower())
459 .set(RunSettingsId::RUN_END_TIME, timeRange.upper())
460 .set(RunSettingsId::RUN_OUTPUT_INTERVAL, timeRange.size() / 10)
461 .set(RunSettingsId::RUN_OUTPUT_TYPE, IoEnum::NONE)
462 .set(RunSettingsId::RUN_OUTPUT_NAME, getIdentifier(name) + "_%d.ssf")
463 .set(RunSettingsId::RUN_VERBOSE_NAME, getIdentifier(name) + ".log")
464 .set(RunSettingsId::SPH_FINDER, FinderEnum::KD_TREE)
465 .set(RunSettingsId::GRAVITY_SOLVER, GravityEnum::BARNES_HUT)
466 .set(RunSettingsId::GRAVITY_KERNEL, GravityKernelEnum::SOLID_SPHERES)
467 .set(RunSettingsId::GRAVITY_OPENING_ANGLE, 0.8_f)
468 .set(RunSettingsId::FINDER_LEAF_SIZE, 20)
469 .set(RunSettingsId::COLLISION_HANDLER, CollisionHandlerEnum::MERGE_OR_BOUNCE)
470 .set(RunSettingsId::COLLISION_OVERLAP, OverlapEnum::PASS_OR_MERGE)
471 .set(RunSettingsId::COLLISION_RESTITUTION_NORMAL, 0.5_f)
472 .set(RunSettingsId::COLLISION_RESTITUTION_TANGENT, 1._f)
473 .set(RunSettingsId::COLLISION_ALLOWED_OVERLAP, 0.01_f)
474 .set(RunSettingsId::COLLISION_BOUNCE_MERGE_LIMIT, 4._f)
475 .set(RunSettingsId::COLLISION_ROTATION_MERGE_LIMIT, 1._f)
476 .set(RunSettingsId::NBODY_INERTIA_TENSOR, false)
477 .set(RunSettingsId::NBODY_MAX_ROTATION_ANGLE, 0.01_f)
478 .set(RunSettingsId::RUN_THREAD_GRANULARITY, 100);
479 return settings;
480 }
481
getSettings()482 VirtualSettings NBodyJob::getSettings() {
483 VirtualSettings connector;
484 addGenericCategory(connector, instName);
485 addTimeSteppingCategory(connector, settings, isResumed);
486 addGravityCategory(connector, settings);
487
488 VirtualSettings::Category& aggregateCat = connector.addCategory("Aggregates (experimental)");
489 aggregateCat.connect<bool>("Enable aggregates", settings, RunSettingsId::NBODY_AGGREGATES_ENABLE);
490 aggregateCat.connect<EnumWrapper>("Initial aggregates", settings, RunSettingsId::NBODY_AGGREGATES_SOURCE)
491 .setEnabler([this] { return settings.get<bool>(RunSettingsId::NBODY_AGGREGATES_ENABLE); });
492
493 VirtualSettings::Category& softCat = connector.addCategory("Soft-body physics (experimental)");
494 softCat.connect("Enable soft-body", "soft.enable", useSoft);
495 softCat.connect<Float>("Repel force strength", settings, RunSettingsId::SOFT_REPEL_STRENGTH)
496 .setEnabler([this] { return useSoft; });
497 softCat.connect<Float>("Friction force strength", settings, RunSettingsId::SOFT_FRICTION_STRENGTH)
498 .setEnabler([this] { return useSoft; });
499
500 auto collisionEnabler = [this] {
501 return !useSoft && !settings.get<bool>(RunSettingsId::NBODY_AGGREGATES_ENABLE) &&
502 settings.get<CollisionHandlerEnum>(RunSettingsId::COLLISION_HANDLER) !=
503 CollisionHandlerEnum::NONE;
504 };
505 auto mergeLimitEnabler = [this] {
506 if (useSoft) {
507 return false;
508 }
509 const CollisionHandlerEnum handler =
510 settings.get<CollisionHandlerEnum>(RunSettingsId::COLLISION_HANDLER);
511 if (handler == CollisionHandlerEnum::NONE) {
512 return false;
513 }
514 const bool aggregates = settings.get<bool>(RunSettingsId::NBODY_AGGREGATES_ENABLE);
515 const OverlapEnum overlap = settings.get<OverlapEnum>(RunSettingsId::COLLISION_OVERLAP);
516 return aggregates || handler == CollisionHandlerEnum::MERGE_OR_BOUNCE ||
517 overlap == OverlapEnum::PASS_OR_MERGE || overlap == OverlapEnum::REPEL_OR_MERGE;
518 };
519
520 VirtualSettings::Category& collisionCat = connector.addCategory("Collisions");
521 collisionCat.connect<EnumWrapper>("Collision handler", settings, RunSettingsId::COLLISION_HANDLER)
522 .setEnabler([this] { //
523 return !useSoft && !settings.get<bool>(RunSettingsId::NBODY_AGGREGATES_ENABLE);
524 });
525 collisionCat.connect<EnumWrapper>("Overlap handler", settings, RunSettingsId::COLLISION_OVERLAP)
526 .setEnabler(collisionEnabler);
527 collisionCat.connect<Float>("Normal restitution", settings, RunSettingsId::COLLISION_RESTITUTION_NORMAL)
528 .setEnabler(collisionEnabler);
529 collisionCat
530 .connect<Float>("Tangential restitution", settings, RunSettingsId::COLLISION_RESTITUTION_TANGENT)
531 .setEnabler(collisionEnabler);
532 collisionCat.connect<Float>("Merge velocity limit", settings, RunSettingsId::COLLISION_BOUNCE_MERGE_LIMIT)
533 .setEnabler(mergeLimitEnabler);
534 collisionCat
535 .connect<Float>("Merge rotation limit", settings, RunSettingsId::COLLISION_ROTATION_MERGE_LIMIT)
536 .setEnabler(mergeLimitEnabler);
537
538 addLoggerCategory(connector, settings);
539 addOutputCategory(connector, settings, *this);
540 return connector;
541 }
542
getRun(const RunSettings & overrides) const543 AutoPtr<IRun> NBodyJob::getRun(const RunSettings& overrides) const {
544 RunSettings run = overrideSettings(settings, overrides, isResumed);
545 return makeAuto<NBodyRun>(run, useSoft);
546 }
547
548 static JobRegistrar sRegisterNBody(
549 "N-body run",
550 "simulations",
__anon0f2951f41e02(const String& name) 551 [](const String& name) { return makeAuto<NBodyJob>(name, EMPTY_SETTINGS); },
552 "Runs N-body simulation using given initial conditions.");
553
554
555 NAMESPACE_SPH_END
556