Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/src/bem_kato.cpp
caa41c47ab18archive/legacy_mpi

Last Change
Raw File

bem_kato.cpp
View Options

/**
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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->functional->computeF();
Real res = this->functional->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->functional->computeGradF();
this->backupTractions();
Real f = 1.;
Real fPrevious = 1e300;
this->nb_iterations = 0;
this->convergence_iterations.clear();
while(f > epsilon && this->nb_iterations < this->max_iterations) {
fPrevious = f;
this->computeDisplacements();
this->functional->computeGradF();
this->backupTractions();
try {
f = linescan(1.,pressure);
}
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
<< 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->functional->getGradF()(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");
UInt n = surface.size();
UInt size = n*n;
//shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (UInt 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