Homec4science
Diffusion tamaas (master) Browse

tamaas/src/bem_kato.cpp
fbfcad1b6cb9master

Last Change
Raw File

bem_kato.cpp
View Options

#include <vector>
#include "surface.hh"
#include "bem_kato.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
Real BemKato::linescan(Real scale,Real pressure){
updatePressure(scale);
shiftPressure(pressure);
truncatePressure();
this->computeF();
Real res = this->getF();
return res;
}
/* -------------------------------------------------------------------------- */
Real BemKato::linesearch(Real & hmax, Real fold,
Real pressure, int search_flag){
if (search_flag){
Real h = hmax;
// Real fold = bem.computeF();
//if (fold == 0) fold = 1e300;
Real f = linescan(h,pressure);
if (f < fold) return f;
while (f > fold){
h *= 0.5;
if (h < 1e-3) throw 1;
f = linescan(h,pressure);
}
f = linescan(h,pressure);
// if (hmax / h > 10) hmax /=2;
return f;
}
return linescan(hmax,pressure);
}
/* -------------------------------------------------------------------------- */
Real BemKato::computeEquilibrium(Real epsilon,
Real pressure,
int nthreads){
this->setNumberOfThreads(nthreads);
// UInt n = surface.size();
// UInt size = n*n;
this->computeDisplacements();
this->computeGradF();
this->backupTractions();
Real f = 1.;
Real fPrevious = 1e300;
Real fOldVer = 1.;
Real fOldVerPrevious = 1e300;
this->nb_iterations = 0;
this->convergence_iterations.clear();
while(f > epsilon && this->nb_iterations < this->max_iterations) {
fPrevious = f;
fOldVerPrevious = fOldVer;
this->computeDisplacements();
this->computeGradF();
this->backupTractions();
try {
f = linescan(1.,pressure);
this->computeOldF();
fOldVer = this->getOldF();
}
catch (int e){
std::cout << " line search failed " << std::endl;
f = linescan(1,pressure);
nb_iterations = 0;
break;
}
if (nb_iterations % dump_freq == 0)
// std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fold << " " << ((f-fold)/forigin) << std::endl;
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " "
<< fOldVer << " " << ((fOldVer-fOldVerPrevious)/fOldVer)
<< std::endl;
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computeTrueDisplacements();
this->computeGaps();
return f;
}
/* -------------------------------------------------------------------------- */
void BemKato::updatePressure(Real scale){
STARTTIMER("updatePressure");
unsigned int n = surface.size();
unsigned int size = n*n;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->surface_tractions(i) =
this->surface_tractions_backup(i) - this->gradient_F(i) * scale;
}
STOPTIMER("updatePressure");
}
/* -------------------------------------------------------------------------- */
void BemKato::backupTractions(){
STARTTIMER("switchPressure");
unsigned int n = surface.size();
unsigned int size = n*n;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->surface_tractions_backup(i) = this->surface_tractions(i);
}
STOPTIMER("switchPressure");
}
/* -------------------------------------------------------------------------- */
Real BemKato::positivePressure(Real step){
STARTTIMER("positivePressure");
unsigned int n = surface.size();
unsigned int size = n*n;
Real p_tot = 0.0;
#pragma omp parallel for reduction(+:p_tot)
for (unsigned int i = 0; i < size; ++i) {
Real sh_press = this->surface_tractions(i) + step;
if (sh_press > max_pressure) p_tot += max_pressure;
else p_tot += sh_press*(sh_press > 0);
}
STOPTIMER("positivePressure");
return p_tot/n/n;
}
/* -------------------------------------------------------------------------- */
void BemKato::shiftPressure(Real required_pressure){
Real step_min = -10;
Real step_max = 10;
STARTTIMER("shiftPressureInitialGuess");
Real p_max = positivePressure(step_max);
Real p_min = positivePressure(step_min);
for(UInt i = 0; p_max < required_pressure && i < 8; ++i, step_max*=10) {
p_max = positivePressure(step_max);
}
for(UInt i = 0; p_min > required_pressure && i < 8; ++i, step_min*=10) {
p_min = positivePressure(step_min);
}
Real p = positivePressure(0.0);
Real epsilon = 1e-12;
STOPTIMER("shiftPressureInitialGuess");
STARTTIMER("shiftPressureDichotomy");
while (fabs(step_min - step_max) > epsilon){
Real step = (step_min + step_max)/2;
p = positivePressure(step);
if (p > required_pressure) step_max = step;
else if (p < required_pressure) step_min = step;
else {
step_max = step;
step_min = step;
}
}
STOPTIMER("shiftPressureDichotomy");
STARTTIMER("shiftPressureTrueShift");
unsigned int n = surface.size();
unsigned int size = n*n;
//shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->surface_tractions(i) += (step_max+step_min)/2;
}
STOPTIMER("shiftPressureTrueShift");
}
/* -------------------------------------------------------------------------- */
void BemKato::truncatePressure(){
STARTTIMER("truncatePressure");
unsigned int n = surface.size();
unsigned int size = n*n;
//shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->surface_tractions(i) > max_pressure)
this->surface_tractions(i) = max_pressure;
else this->surface_tractions(i) *= (this->surface_tractions(i) > 0);
}
STOPTIMER("truncatePressure");
}
/* -------------------------------------------------------------------------- */
c4science · Help