Homec4science
Diffusion yearII_reports (master) Browse

rSTICAZZI/physnum/rap4/Exercice3.cpp
f2a8d046e9ebmaster

Last Change
Raw File

Exercice3.cpp
View Options

#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
#include "ConfigFile.tpp" // Fichier .tpp car inclut un template
using namespace std;
// verlet schema implementation
// x = position
// v = speed
// t = time
// a = acceleration function, owned by a Base class
// returns the delta of position dx
template<typename T, class Base>
void verlet_integration(T &x, T &v, const double t, const double dt, T (Base::*a)(T, T, double) const, Base *base)
{
T x_old(x), v_old(v);
#define ACC(x, v, t) (base->*a)(x, v, t)
x += v * dt + 0.5 * ACC(x, v, t) * dt*dt;
v += 0.5 * ACC(x_old, v_old ,t) *dt;
v = v_old + 0.5 * dt * ( ACC(x_old, v, t) + ACC(x, v, t+dt) );
}
double bound(double theta)
{
while(theta > M_PI || theta < - M_PI)
theta -= 2 * M_PI * ((theta > M_PI) ? 1 : -1);
return theta;
}
struct Errorvaluation
{
Vector v;
double err;
};
class Sym
{
private:
#define CONVN 4
double t, dt, tFin;
double rho_0, S, C_x, Rt, lambda, G, eps;
int sampling;
int last;
size_t count;
ofstream *outputFile;
void printOut(bool force)
{
if((!force && last>=sampling) || (force && last!=1))
{
double emec = kinetic_energy() + U();
if (!exo3eII || (count >= countoff * N_e && (count % N_e) == 0)) // only for exercise 3.3e)
*outputFile << t << " " << bound(theta) << " " << thetadot << " " << emec << " " << pnc() << " " << pc() << " " << pc() - pnc() << endl;
last = 1;
}
else
{
last++;
}
}
double kinetic_energy() const
{
double v_r = L * thetadot;
return 0.5 * m * v_r * v_r;
}
double U() const
{
double h_r = L * (1 - cos(theta));
return m * g * h_r;
}
double pnc() const
{
return
- kappa*L*L*thetadot*thetadot
+ m*L*Omega*Omega*d* sin(Omega*t) * sin(theta) * thetadot;
}
double pc() const
{
return
m * L* L * thetadot * acctet(theta, thetadot, t)
+ m * g * L * sin(theta) * thetadot;
}
double acctet(double thet, double thetad, double temps) const {
double a(g*sin(thet)/L), b(kappa*thetad/m), c(Omega*Omega*d/L), e(sin(Omega*temps)*sin(thet));
return -a-b+c*e;
}
void step()
{
verlet_integration<double, Exercice3>(theta, thetadot, t, dt, &Exercice3::acctet, this);
count++;
}
/*
Errorvaluation error(const Vector & argum, double dt, Vector(Sym::*evolution)(const Vector&, double), Sym *sym){
Vector y1 = runge_kutta<Vector, Sym>(argum, dt, 0, evolution, sym);
Vector ys = runge_kutta<Vector, Sym>(argum, dt/2, 0, evolution, sym);
Vector y2 = runge_kutta<Vector, Sym>(ys, dt/2, 0, evolution, sym);
Vector yf = y1 - y2;
Errorvaluation e;
e.err = yf.sq_module();
e.v =y2;
return e;
}
double changedt(double dt,double d){
return 0.99*dt*pow(eps/d, 1/(CONVN+1));
}
void adaptifstep(Vector &argum, double &dt, Vector(Sym::*evolution)(const Vector&, double), Sym *sym){
Errorvaluation e = error(argum, dt, evolution, sym);
double d = e.err;
while (d >= eps){
dt = changedt(dt,d);
e = error(argum, dt, evolution, sym);
d = e.err;
}
argum=e.v;
dt = changedt(dt,d);
}
*/
public:
Sym(int argc, char* argv[])
{
string inputPath("configuration.in"); // Fichier d'input par defaut
if(argc>1) // Fichier d'input specifie par l'utilisateur ("./Exercice3 config_perso.in")
inputPath = argv[1];
ConfigFile configFile(inputPath); // Les parametres sont lus et stockes dans une "map" de strings.
for(int i(2); i<argc; ++i) // Input complementaires ("./Exercice3 config_perso.in input_scan=[valeur]")
configFile.process(argv[i]);
tFin = configFile.get<double>("tFin");
dt = configFile.get<double>("dt");
d = configFile.get<double>("d");
Omega = configFile.get<double>("Omega");
kappa = configFile.get<double>("kappa");
m = configFile.get<double>("m");
g = configFile.get<double>("g");
L = configFile.get<double>("L");
theta = configFile.get<double>("theta0");
thetadot = configFile.get<double>("thetadot0");
sampling = configFile.get<int>("sampling");
exo3eII = configFile.get<int>("exo3eII") == 1;
countoff = configFile.get<int>("countoff");
// Ouverture du fichier de sortie
outputFile = new ofstream(configFile.get<string>("output").c_str());
outputFile->precision(15);
count = 0;
N_e = (2 * M_PI / (Omega * dt));
};
~Exercice3()
{
outputFile->close();
delete outputFile;
};
void run()
{
t = 0.;
last = 0;
printOut(true);
while( t < tFin-0.5*dt )
{
step();
t += dt;
printOut(false);
}
printOut(true);
};
};
int main(int argc, char* argv[])
{
Sym engine(argc, argv);
engine.run();
return 0;
}
c4science · Help