bubble-dynamics

bubbleDynamics.cpp (6364B)

---

 /* File:   bubbleDynamics.cpp */
 /* Date:   Tue Jan 12 22:04:58 2016 */
 /* Author: Fabian Wermelinger */
 /* Tag:    Solve 1D bubble dynamics */
 /* Copyright 2016 ETH Zurich. All Rights Reserved. */
 #include "kernels.h"
 
 #include <ODETB/GnuplotDump.h>
 #include <ODETB/TimeStepper/TimeStepper.h>
 #include <ODETB/common.h>
 
 #include <cassert>
 #include <cmath>
 #include <cstdio>
 #include <ctime>
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <memory>
 #include <stdexcept>
 
 using namespace std;
 
 
 void say_hello(const int argc, const char * argv[])
 {
     time_t rawtime;
     struct tm* timeinfo;
     time(&rawtime);
     timeinfo = localtime(&rawtime);
     char buf[256];
     strftime(buf, 256, "%A, %h %d %Y, %r", timeinfo);
     cout << "Hi! It is now " << buf << endl;
     cout << "Here is what you gave me: ";
     for (int i = 0; i < argc; ++i)
         cout << argv[i] << " ";
     cout << endl;
 }
 
 template <typename T>
 void run(ArgumentParser &parser,
          StepperSettings &stepperData,
          BubbleData &simConfig,
          StepperBase<T> *const stepper)
 {
     simConfig.print();
 
     using vec_t = T;
     T &U = stepper->U();
 
     size_t& step       = stepperData.step;
     Real& t            = stepperData.t;
     Real& dt           = stepperData.dt;
     const Real& tFinal = stepperData.tFinal;
     Real& tDump        = stepperData.tDump;
     if (stepperData.dtDump > 0.0) tDump = t + stepperData.dtDump;
 
 
     if (stepperData.bRestart)
     {
         deserialize<vec_t>(U, stepperData);
         std::cout << "-->RESTARING: Time = " << std::scientific << stepperData.t;
         std::cout << ";\tStep = " << std::fixed << stepperData.step << std::endl;
     }
     else
         stepper->write(step, t, dt, U, &simConfig);
 
 
     while (t < tFinal)
     {
         dt = (tFinal-t) < dt ? (tFinal-t) : dt;
         if (stepperData.dtDump > 0.0) dt = (tDump-t) < dt ? (tDump-t) : dt;
 
         stepper->step(&simConfig);
 
         if (stepperData.bFixedStep && stepperData.dtDump > 0.0)
         {
             if (t == tDump)
             {
                 stepper->write(step, t, dt, U, &simConfig);
                 tDump += stepperData.dtDump;
             }
         }
         else
             if (stepperData.bFixedStep && step % stepperData.writeGranularity == 0)
                 stepper->write(step, t, dt, U, &simConfig);
 
         if (stepperData.restartstep > 0 && step % stepperData.restartstep == 0)
             serialize<vec_t>(U, stepperData);
     }
 
     if (stepperData.dtDump > 0.0)
     {
         if (t != tDump)
             stepper->write(step, t, dt, U, &simConfig);
     }
     else
     {
         if (step % stepperData.writeGranularity != 0)
             stepper->write(step, t, dt, U, &simConfig);
     }
 
     printf("\nFinal Time = %e;\nFinal Step = %zu\n",
            stepperData.t,
            stepperData.step);
 }
 
 int main(int argc, const char **argv)
 {
     say_hello(argc, argv);
 
     // Parameter
     ArgumentParser parser(argc, argv);
     if (parser.exist("-conf"))
         parser.readFile(parser("-conf").asString());
     parser.print_args();
 
     // set user specific data
     StepperSettings stepperData(parser);
     BubbleData simConfig(parser);
 
     int binary_dump;
     if (parser("-binary").asBool(false))
         binary_dump = GnuplotDump::BIN;
     else
         binary_dump = GnuplotDump::ASCII;
 
     if (parser("-dynamics").asString("rp") == "rpcluster-positions" ||
         parser("-dynamics").asString("rp") == "kmcluster-positions") {
         // with additional 6 DoF for bubble translation (3 x and 3 \dot{x})
         using TVec = LightVector<State<Real, 8>>;
 
         // Select stepper and kernel
         TVec U(simConfig.Nbubbles);
         std::unique_ptr<KernelBase<TVec>> kern;
         if (parser("-dynamics").asString("rp") == "rpcluster-positions") {
             RPClusterPositions_D<TVec>::IC(U, simConfig);
             kern.reset(new RPClusterPositions_D<TVec>(U.size(), binary_dump));
         } else if (parser("-dynamics").asString("rp") ==
                    "kmcluster-positions") {
             KMClusterPositions_D<TVec>::IC(U, simConfig);
             kern.reset(new KMClusterPositions_D<TVec>(U.size(), binary_dump));
         } else {
             throw std::runtime_error("Unknown kernel: " +
                                      parser("-dynamics").asString("rp"));
         }
 
         std::unique_ptr<StepperBase<TVec>> stepper(
             stepperData.template stepperFactory<TVec>(parser, U, kern.get()));
         assert(stepper != nullptr);
 
         run<TVec>(parser, stepperData, simConfig, stepper.get());
     } else {
         // regular R and \dot{R} state
         using TVec = LightVector<State<Real, 2>>;
 
         // Select stepper and kernel
         TVec U(simConfig.Nbubbles);
         std::unique_ptr<KernelBase<TVec>> kern;
         if (parser("-dynamics").asString("rp") == "rp") {
             RayleighPlesset<TVec>::IC(U, simConfig);
             kern.reset(new RayleighPlesset<TVec>(U.size(), binary_dump));
         } else if (parser("-dynamics").asString("rp") == "hbgl") {
             RayleighPlesset_HBGL<TVec>::IC(U, simConfig);
             kern.reset(new RayleighPlesset_HBGL<TVec>(U.size(), binary_dump));
         } else if (parser("-dynamics").asString("rp") == "km") {
             KellerMiksis<TVec>::IC(U, simConfig);
             kern.reset(new KellerMiksis<TVec>(U.size(), binary_dump));
         } else if (parser("-dynamics").asString("rp") == "kmcluster") {
             KMCluster_TY<TVec>::IC(U, simConfig);
             kern.reset(new KMCluster_TY<TVec>(U.size(), binary_dump));
         } else if (parser("-dynamics").asString("rp") == "kmcluster-fc") {
             KMCluster_FC<TVec>::IC(U, simConfig);
             kern.reset(new KMCluster_FC<TVec>(U.size(), binary_dump));
         } else if (parser("-dynamics").asString("rp") == "rpcluster") {
             RPCluster<TVec>::IC(U, simConfig);
             kern.reset(new RPCluster<TVec>(U.size(), binary_dump));
         } else {
             throw std::runtime_error("Unknown kernel: " +
                                      parser("-dynamics").asString("rp"));
         }
 
         std::unique_ptr<StepperBase<TVec>> stepper(
             stepperData.template stepperFactory<TVec>(parser, U, kern.get()));
         assert(stepper != nullptr);
 
         run<TVec>(parser, stepperData, simConfig, stepper.get());
     }
     return 0;
 }