1 /// \brief Executable running a single impact simulation, using command-line parameters
2
3 #include "Sph.h"
4 #include "run/Node.h"
5 #include "run/jobs/InitialConditionJobs.h"
6 #include "run/jobs/IoJobs.h"
7 #include "run/jobs/MaterialJobs.h"
8 #include "run/jobs/ParticleJobs.h"
9 #include "run/jobs/SimulationJobs.h"
10
11 using namespace Sph;
12
13 static Array<ArgDesc> params{
14 { "tr", "target-radius", ArgEnum::FLOAT, "Radius of the target [m]" },
15 { "tp", "target-period", ArgEnum::FLOAT, "Rotational period of the target [h]" },
16 { "ir", "impactor-radius", ArgEnum::FLOAT, "Radius of the impactor [m]" },
17 { "q",
18 "impact-energy",
19 ArgEnum::FLOAT,
20 "Relative impact energy Q/Q_D^*. This option can be only used together with -tr and -v, it is "
21 "incompatible with -ir." },
22 { "v", "impact-speed", ArgEnum::FLOAT, "Impact speed [km/s]" },
23 { "phi", "impact-angle", ArgEnum::FLOAT, "Impact angle [deg]" },
24 { "n", "particle-count", ArgEnum::INT, "Number of particles in the target" },
25 { "st", "stabilization-time", ArgEnum::FLOAT, "Duration of the stabilization phase [s]" },
26 { "ft", "fragmentation-time", ArgEnum::FLOAT, "Duration of the fragmentation phase [s]" },
27 { "rt", "reaccumulation-time", ArgEnum::FLOAT, "Duration of the reaccumulation phase [s]" },
28 { "i", "resume-from", ArgEnum::STRING, "Resume simulation from given state file" },
29 { "o",
30 "output-dir",
31 ArgEnum::STRING,
32 "Directory containing configuration files and run output files. If not specified, it is determined "
33 "from other arguments. If no arguments are specified, the current working directory is used." },
34 };
35
printBanner(ILogger & logger)36 static void printBanner(ILogger& logger) {
37 logger.write("*******************************************************************************");
38 logger.write("******************************* OpenSPH Impact ********************************");
39 logger.write("*******************************************************************************");
40 }
41
printNoConfigsMsg(ILogger & logger,const Path & outputDir)42 static void printNoConfigsMsg(ILogger& logger, const Path& outputDir) {
43 logger.write("");
44 logger.write("No configuration files found, the program will generate default configuration ");
45 logger.write("files and save them to directory '" + outputDir.string() + "'");
46 logger.write("");
47 logger.write("To start a simulation, re-run this program; it will load the generated files. ");
48 logger.write("You can also specify parameters of the simulation as command-line arguments. ");
49 logger.write("Note that these arguments will override parameters loaded from configuration ");
50 logger.write("files. For more information, execute the program with -h (or --help) argument.");
51 logger.write("");
52 }
53
54 template <typename TEnum>
loadSettings(const Path & path,const Settings<TEnum> & defaults,ILogger & logger,bool & doDryRun)55 static Settings<TEnum> loadSettings(const Path& path,
56 const Settings<TEnum>& defaults,
57 ILogger& logger,
58 bool& doDryRun) {
59 Settings<TEnum> settings = defaults;
60 const bool result = settings.tryLoadFileOrSaveCurrent(path);
61 if (result) {
62 logger.write("Loaded configuration file '" + path.string() + "'");
63
64 // at least one configuration file exists, run the simulation
65 doDryRun = false;
66 } else {
67 logger.write("No file '" + path.string() + "' found, it has been created with default parameters");
68 }
69 return settings;
70 }
71
72 struct RunParams {
73 Optional<Float> targetRadius;
74 Optional<Float> targetPeriod;
75 Optional<Float> impactorRadius;
76 Optional<Float> impactAngle;
77 Optional<Float> impactSpeed;
78 Optional<int> particleCnt;
79
80 Optional<Float> stabTime;
81 Optional<Float> fragTime;
82 Optional<Float> reacTime;
83
84 Optional<String> resumePath;
85 Optional<String> outputPath;
86
87
88 /// \brief Returns the name of the created output directory.
getOutputPathRunParams89 String getOutputPath() const {
90 if (outputPath) {
91 return outputPath.value();
92 }
93
94 std::stringstream ss;
95 ss << "sph_";
96 if (targetRadius) {
97 ss << round(targetRadius.value()) << "m_";
98 }
99 if (impactorRadius) {
100 ss << round(impactorRadius.value()) << "m_";
101 }
102 if (targetPeriod) {
103 ss << round(60._f * targetPeriod.value()) << "min_";
104 }
105 if (impactSpeed) {
106 ss << round(impactSpeed.value() / 1.e3_f) << "kms_";
107 }
108 if (impactAngle) {
109 ss << round(impactAngle.value()) << "deg_";
110 }
111 if (particleCnt) {
112 ss << particleCnt.value() << "p_";
113 }
114
115 String name = String::fromAscii(ss.str().c_str());
116 // drop the last "_";
117 name.erase(name.size() - 1, 1);
118 return name;
119 }
120 };
121
122 class RunFactory {
123 private:
124 StringLogger logger;
125 RunParams params;
126 Path outputDir;
127 bool doDryRun;
128 String paramsMsg;
129
130 public:
RunFactory(const RunParams & params)131 RunFactory(const RunParams& params)
132 : params(params) {
133 doDryRun = true;
134 outputDir = Path(params.getOutputPath());
135 }
136
makeSimulation()137 SharedPtr<JobNode> makeSimulation() {
138 if (params.resumePath) {
139 BinaryInput input;
140 Expected<BinaryInput::Info> info = input.getInfo(Path(params.resumePath.value()));
141 if (!info) {
142 throw Exception("Cannot resume simulation from file '" + params.resumePath.value() + "'.\n" +
143 info.error());
144 }
145
146 logger.write("Resuming simulation from file '" + params.resumePath.value() + "'");
147 if (info->runType == RunTypeEnum::SPH) {
148 return resumeFragmentation();
149 } else if (info->runType == RunTypeEnum::NBODY) {
150 return resumeReaccumulation();
151 } else {
152 throw Exception("Cannot resume from this file.");
153 }
154
155 } else {
156 logger.write("Starting new simulation");
157 return makeNewSimulation();
158 }
159 }
160
isDryRun() const161 bool isDryRun() const {
162 return doDryRun;
163 }
164
getBannerMsg() const165 String getBannerMsg() const {
166 return logger.toString() + "\n" + paramsMsg;
167 }
168
getOutputDir() const169 Path getOutputDir() const {
170 return outputDir;
171 }
172
173 private:
defaultSphTime(const Optional<Float> runTime,const Optional<Float> radius,const Float mult)174 Float defaultSphTime(const Optional<Float> runTime, const Optional<Float> radius, const Float mult) {
175 // by default, use 1h for 100km body
176 return runTime.valueOr(mult * 3600._f * radius.valueOr(5.e4_f) / 5.e4_f);
177 }
178
overrideRunTime(RunSettings & settings,const Float endTime)179 void overrideRunTime(RunSettings& settings, const Float endTime) {
180 settings.set(RunSettingsId::RUN_END_TIME, endTime)
181 .set(RunSettingsId::RUN_OUTPUT_INTERVAL, endTime / 10._f);
182 }
183
makeCollisionSetup()184 SharedPtr<JobNode> makeCollisionSetup() {
185 // target IC
186 BodySettings targetDefaults;
187 targetDefaults.set(BodySettingsId::BODY_RADIUS, params.targetRadius.valueOr(50.e3_f))
188 .set(BodySettingsId::PARTICLE_COUNT, params.particleCnt.valueOr(10000));
189 if (params.targetPeriod) {
190 targetDefaults.set(BodySettingsId::BODY_SPIN_RATE, 24._f / params.targetPeriod.value());
191 }
192 BodySettings targetBody =
193 loadSettings<BodySettingsId>(outputDir / Path("target.cnf"), targetDefaults, logger, doDryRun);
194 SharedPtr<JobNode> targetIc = makeNode<MonolithicBodyIc>("target body", targetBody);
195
196 // impactor IC
197 BodySettings impactorDefaults;
198 impactorDefaults.set(BodySettingsId::BODY_RADIUS, params.targetRadius.valueOr(10.e3_f))
199 .set(BodySettingsId::DAMAGE_MIN, LARGE)
200 .set(BodySettingsId::STRESS_TENSOR_MIN, LARGE);
201 impactorDefaults.unset(BodySettingsId::PARTICLE_COUNT); // never used
202 BodySettings impactorBody = loadSettings<BodySettingsId>(
203 outputDir / Path("impactor.cnf"), impactorDefaults, logger, doDryRun);
204 SharedPtr<JobNode> impactorIc = makeNode<ImpactorIc>("impactor body", impactorBody);
205 targetIc->connect(impactorIc, "target");
206
207 // target stabilization
208 RunSettings stabDefaults = SphStabilizationJob::getDefaultSettings("stabilization");
209 stabDefaults.set(RunSettingsId::RUN_OUTPUT_PATH, outputDir.string());
210 overrideRunTime(stabDefaults, defaultSphTime(params.stabTime, params.targetRadius, 0.2_f));
211
212 RunSettings stabRun =
213 loadSettings<RunSettingsId>(outputDir / Path("stab.cnf"), stabDefaults, logger, doDryRun);
214 SharedPtr<JobNode> stabTarget = makeNode<SphStabilizationJob>("stabilization", stabRun);
215 targetIc->connect(stabTarget, "particles");
216
217 // collision setup
218 CollisionGeometrySettings geometryDefaults;
219 geometryDefaults.set(CollisionGeometrySettingsId::IMPACT_SPEED, params.impactSpeed.valueOr(5.e3_f))
220 .set(CollisionGeometrySettingsId::IMPACT_ANGLE, params.impactAngle.valueOr(45._f));
221 CollisionGeometrySettings geometry = loadSettings<CollisionGeometrySettingsId>(
222 outputDir / Path("geometry.cnf"), geometryDefaults, logger, doDryRun);
223 SharedPtr<JobNode> setup = makeNode<CollisionGeometrySetup>("geometry", geometry);
224 stabTarget->connect(setup, "target");
225 impactorIc->connect(setup, "impactor");
226
227 printRunSettings(targetBody, impactorBody, geometry);
228
229 return setup;
230 }
231
makeFragmentation()232 SharedPtr<JobNode> makeFragmentation() {
233 RunSettings fragDefaults = SphJob::getDefaultSettings("fragmentation");
234 fragDefaults.set(RunSettingsId::RUN_OUTPUT_PATH, outputDir.string());
235 overrideRunTime(fragDefaults, defaultSphTime(params.fragTime, params.targetRadius, 1._f));
236 RunSettings fragRun =
237 loadSettings<RunSettingsId>(outputDir / Path("frag.cnf"), fragDefaults, logger, doDryRun);
238 SharedPtr<JobNode> frag = makeNode<SphJob>("fragmentation", fragRun);
239 return frag;
240 }
241
makeReaccumulation()242 SharedPtr<JobNode> makeReaccumulation() {
243 RunSettings reacDefaults = NBodyJob::getDefaultSettings("reaccumulation");
244 reacDefaults.set(RunSettingsId::RUN_OUTPUT_PATH, outputDir.string());
245 overrideRunTime(reacDefaults, params.reacTime.valueOr(3600._f * 24._f * 10._f));
246 RunSettings reacRun =
247 loadSettings<RunSettingsId>(outputDir / Path("reac.cnf"), reacDefaults, logger, doDryRun);
248 SharedPtr<JobNode> reac = makeNode<NBodyJob>("reaccumulation", reacRun);
249 return reac;
250 }
251
makeNewSimulation()252 SharedPtr<JobNode> makeNewSimulation() {
253 SharedPtr<JobNode> setup = makeCollisionSetup();
254 SharedPtr<JobNode> frag = makeFragmentation();
255 setup->connect(frag, "particles");
256
257 SharedPtr<JobNode> handoff = makeNode<SmoothedToSolidHandoff>("handoff"); // has no parameters
258 frag->connect(handoff, "particles");
259
260 SharedPtr<JobNode> reac = makeReaccumulation();
261 handoff->connect(reac, "particles");
262 return reac;
263 }
264
resumeFragmentation()265 SharedPtr<JobNode> resumeFragmentation() {
266 SharedPtr<JobNode> loadFile = makeNode<LoadFileJob>(Path(params.resumePath.value()));
267
268 SharedPtr<JobNode> frag = makeFragmentation();
269 loadFile->connect(frag, "particles");
270
271 SharedPtr<JobNode> handoff = makeNode<SmoothedToSolidHandoff>("handoff");
272 frag->connect(handoff, "particles");
273
274 SharedPtr<JobNode> reac = makeReaccumulation();
275 handoff->connect(reac, "particles");
276 return reac;
277 }
278
resumeReaccumulation()279 SharedPtr<JobNode> resumeReaccumulation() {
280 SharedPtr<JobNode> loadFile = makeNode<LoadFileJob>(Path(params.resumePath.value()));
281
282 SharedPtr<JobNode> reac = makeReaccumulation();
283 loadFile->connect(reac, "particles");
284 return reac;
285 }
286
printRunSettings(const BodySettings & targetBody,const BodySettings & impactorBody,const CollisionGeometrySettings & geometry)287 void printRunSettings(const BodySettings& targetBody,
288 const BodySettings& impactorBody,
289 const CollisionGeometrySettings& geometry) {
290 const Float targetRadius = targetBody.get<Float>(BodySettingsId::BODY_RADIUS);
291 const Float impactorRadius = impactorBody.get<Float>(BodySettingsId::BODY_RADIUS);
292 const Float impactSpeed = geometry.get<Float>(CollisionGeometrySettingsId::IMPACT_SPEED);
293 const Float impactAngle = geometry.get<Float>(CollisionGeometrySettingsId::IMPACT_ANGLE);
294 const Float spinRate = targetBody.get<Float>(BodySettingsId::BODY_SPIN_RATE);
295 const Size particleCnt = targetBody.get<int>(BodySettingsId::PARTICLE_COUNT);
296 const Float rho = targetBody.get<Float>(BodySettingsId::DENSITY);
297 const Float Q_D = evalBenzAsphaugScalingLaw(2._f * targetRadius, rho);
298 const Float impactEnergy = getImpactEnergy(targetRadius, impactorRadius, impactSpeed) / Q_D;
299
300 StringLogger logger;
301 logger.setScientific(false);
302 logger.setPrecision(4);
303 logger.write();
304 logger.write("Run parameters");
305 logger.write("-------------------------------------");
306 logger.write(" Target radius (R_pb): ", 1.e-3_f * targetRadius, " km");
307 logger.write(" Impactor radius (r_imp): ", 1.e-3_f * impactorRadius, " km");
308 logger.write(" Impact speed (v_imp): ", 1.e-3_f * impactSpeed, " km/s");
309 logger.write(" Impact angle (phi_imp): ", impactAngle, "°");
310 logger.write(" Impact energy (Q/Q*_D): ", impactEnergy);
311 logger.write(" Target period (P_pb): ", 24._f / spinRate, spinRate == 0._f ? "" : "h");
312 logger.write(" Particle count (N): ", particleCnt);
313 logger.write("-------------------------------------");
314 logger.write();
315 logger.setScientific(true);
316 logger.setPrecision(PRECISION);
317
318 paramsMsg = logger.toString();
319 }
320 };
321
run(const ArgParser & parser,ILogger & logger)322 static void run(const ArgParser& parser, ILogger& logger) {
323 RunParams params;
324 params.targetRadius = parser.tryGetArg<Float>("tr");
325 params.targetPeriod = parser.tryGetArg<Float>("tp");
326 params.impactSpeed = parser.tryGetArg<Float>("v");
327 if (params.impactSpeed) {
328 params.impactSpeed.value() *= 1.e3_f; // km/s -> m/s
329 }
330 params.impactAngle = parser.tryGetArg<Float>("phi");
331 params.impactorRadius = parser.tryGetArg<Float>("ir");
332 params.particleCnt = parser.tryGetArg<int>("n");
333
334 if (Optional<Float> impactEnergy = parser.tryGetArg<Float>("q")) {
335 // we have to specify also -tr and -v, as output directory is determined fom the computed
336 // impactorRadius. We cannot use values loaded from config files, as it would create circular
337 // dependency: we need impactor radius to get output path, we need output path to load config
338 // files, we need config files to get impactor radius ...
339 const Float density = BodySettings::getDefaults().get<Float>(BodySettingsId::DENSITY);
340 const Optional<Float> targetRadius = parser.tryGetArg<Float>("tr");
341 const Optional<Float> impactSpeed = parser.tryGetArg<Float>("v");
342 if (!targetRadius || !impactSpeed) {
343 throw ArgError(
344 "To specify impact energy (-q), you also need to specify the target radius (-tr) and "
345 "impact speed (-v)");
346 }
347 params.impactorRadius = getImpactorRadius(
348 targetRadius.value(), impactSpeed.value() * 1.e3_f, impactEnergy.value(), density);
349 }
350
351 params.outputPath = parser.tryGetArg<String>("o");
352 params.resumePath = parser.tryGetArg<String>("i");
353
354 params.stabTime = parser.tryGetArg<Float>("st");
355 params.fragTime = parser.tryGetArg<Float>("ft");
356 params.reacTime = parser.tryGetArg<Float>("rt");
357
358 RunFactory factory(params);
359 SharedPtr<JobNode> node = factory.makeSimulation();
360
361 printBanner(logger);
362 if (!factory.isDryRun()) {
363 logger.write(factory.getBannerMsg());
364
365 NullJobCallbacks callbacks;
366 node->run(EMPTY_SETTINGS, callbacks);
367 } else {
368 printNoConfigsMsg(logger, factory.getOutputDir());
369 }
370 }
371
main(int argc,char * argv[])372 int main(int argc, char* argv[]) {
373 StdOutLogger logger;
374
375 try {
376 ArgParser parser(params);
377 parser.parse(argc, argv);
378
379 run(parser, logger);
380
381 return 0;
382
383 } catch (HelpException& e) {
384 printBanner(logger);
385 logger.write(e.what());
386 return 0;
387 } catch (const Exception& e) {
388 logger.write("Run failed!\n", e.what());
389 return -1;
390 }
391 }
392