Newer
Older
#include <dune/common/bitfield.hh>
#include <dune/istl/io.hh>
#include <dune/disc/functions/p1function.hh>
#include <dune/disc/miscoperators/laplace.hh>
#include <dune/disc/operators/p1operator.hh>
#include "../common/trustregiongsstep.hh"
#include "../contact/src/contactmmgstep.hh"
#include "../solver/iterativesolver.hh"
#include "../common/energynorm.hh"
#include "../common/h1seminorm.hh"
#include "configuration.hh"
#include "quaternion.hh"
// for debugging
#include "../test/fdcheck.hh"
// Number of degrees of freedom:
// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
//const int blocksize = 6;
template <class GridType>
void RodSolver<GridType>::
setTrustRegionObstacles(double trustRegionRadius,
std::vector<BoxConstraint<blocksize> >& trustRegionObstacles)
{
for (int j=0; j<trustRegionObstacles.size(); j++) {
for (int k=0; k<blocksize; k++) {
trustRegionObstacles[j].val[2*k] = -trustRegionRadius;
trustRegionObstacles[j].val[2*k+1] = trustRegionRadius;
}
}
//std::cout << trustRegionObstacles << std::endl;
// exit(0);
}
template <class GridType>
void RodSolver<GridType>::setup(const GridType& grid,
const RodAssembler<GridType>* rodAssembler,
int maxTrustRegionSteps,
double initialTrustRegionRadius,
int multigridIterations,
double mgTolerance,
int mu,
int nu1,
int nu2,
int baseIterations,
double baseTolerance,
bool instrumented)
{
using namespace Dune;
grid_ = &grid;
rodAssembler_ = rodAssembler;
x_ = x;
maxTrustRegionSteps_ = maxTrustRegionSteps;
initialTrustRegionRadius_ = initialTrustRegionRadius;
multigridIterations_ = multigridIterations;
qpTolerance_ = mgTolerance;
mu_ = mu;
nu1_ = nu1;
nu2_ = nu2;
baseIt_ = baseIterations;
baseTolerance_ = baseTolerance;
instrumented_ = instrumented;
int numLevels = grid_->maxLevel()+1;
// ////////////////////////////////////////////
// Construct array with the Dirichlet nodes
// ////////////////////////////////////////////
dirichletNodes_.resize(numLevels);
for (int i=0; i<numLevels; i++) {
dirichletNodes_[i].resize( blocksize * grid.size(i,1), false );
for (int j=0; j<blocksize; j++) {
dirichletNodes_[i][j] = true;
dirichletNodes_[i][dirichletNodes_[i].size()-1-j] = true;
}
}
// ////////////////////////////////
// Create a multigrid solver
// ////////////////////////////////
// First create a gauss-seidel base solver
TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep);
IterativeSolver<MatrixType, CorrectionType>* baseSolver
= new IterativeSolver<MatrixType, CorrectionType>(baseSolverStep,
baseIt_,
baseTolerance_,
baseEnergyNorm,
Solver::QUIET);
// Make pre and postsmoothers
TrustRegionGSStep<MatrixType, CorrectionType>* presmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
TrustRegionGSStep<MatrixType, CorrectionType>* postsmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
ContactMMGStep<MatrixType, CorrectionType>* mmgStep = new ContactMMGStep<MatrixType, CorrectionType>(numLevels);
mmgStep->setMGType(mu_, nu1_, nu2_);
mmgStep->dirichletNodes_ = &dirichletNodes_;
mmgStep->basesolver_ = baseSolver;
mmgStep->presmoother_ = presmoother;
mmgStep->postsmoother_ = postsmoother;
mmgStep->hasObstacle_ = &hasObstacle_;
mmgStep->obstacles_ = &trustRegionObstacles_;
// //////////////////////////////////////////////////////////////////////////////////////
// Assemble a Laplace matrix to create a norm that's equivalent to the H1-norm
// //////////////////////////////////////////////////////////////////////////////////////
LeafP1Function<GridType,double> u(grid),f(grid);
LaplaceLocalStiffness<GridType,double> laplaceStiffness;
LeafP1OperatorAssembler<GridType,double,1>* A = new LeafP1OperatorAssembler<GridType,double,1>(grid);
A->assemble(laplaceStiffness,u,f);
typedef typename LeafP1OperatorAssembler<GridType,double,1>::RepresentationType LaplaceMatrixType;
if (h1SemiNorm_)
delete h1SemiNorm_;
h1SemiNorm_ = new H1SemiNorm<CorrectionType>(**A);
mmgSolver_ = new IterativeSolver<MatrixType, CorrectionType>(mmgStep,
multigridIterations_,
qpTolerance_,
h1SemiNorm_,
// Write all intermediate solutions, if requested
if (instrumented_)
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
// ////////////////////////////////////////////////////////////
// Create Hessian matrix and its occupation structure
// ////////////////////////////////////////////////////////////
if (hessianMatrix_)
delete hessianMatrix_;
hessianMatrix_ = new MatrixType;
MatrixIndexSet indices(grid_->size(1), grid_->size(1));
rodAssembler_->getNeighborsPerVertex(indices);
indices.exportIdx(*hessianMatrix_);
// //////////////////////////////////////////////////////////
// Create obstacles
// //////////////////////////////////////////////////////////
hasObstacle_.resize(numLevels);
for (int i=0; i<hasObstacle_.size(); i++)
hasObstacle_[i].resize(grid_->size(i, 1),true);
trustRegionObstacles_.resize(numLevels);
for (int i=0; i<numLevels; i++)
trustRegionObstacles_[i].resize(grid_->size(i,1));
// ////////////////////////////////////
// Create the transfer operators
// ////////////////////////////////////
for (int k=0; k<mmgStep->mgTransfer_.size(); k++)
delete(mmgStep->mgTransfer_[k]);
mmgStep->mgTransfer_.resize(numLevels-1);
for (int i=0; i<mmgStep->mgTransfer_.size(); i++){
TruncatedMGTransfer<CorrectionType>* newTransferOp = new TruncatedMGTransfer<CorrectionType>;
newTransferOp->setup(*grid_,i,i+1);
mmgStep->mgTransfer_[i] = newTransferOp;
}
}
template <class GridType>
void RodSolver<GridType>::solve()
{
using namespace Dune;
double trustRegionRadius = initialTrustRegionRadius_;
// /////////////////////////////////////////////////////
// Set up the log file, if requested
// /////////////////////////////////////////////////////
FILE* fp;
if (instrumented_) {
fp = fopen("statistics", "w");
if (!fp)
DUNE_THROW(IOError, "Couldn't open statistics file for writing!");
// /////////////////////////////////////////////////////
// Trust-Region Solver
// /////////////////////////////////////////////////////
for (int i=0; i<maxTrustRegionSteps_; i++) {
if (this->verbosity_ == FULL) {
std::cout << "----------------------------------------------------" << std::endl;
std::cout << " Trust-Region Step Number: " << i << std::endl;
std::cout << "----------------------------------------------------" << std::endl;
std::cout << "### Trust-Region Radius: " << trustRegionRadius << " ###" << std::endl;
}
CorrectionType rhs;
CorrectionType corr(x_.size());
corr = 0;
if (this->verbosity_ == FULL)
std::cout << "Rod energy: " <<rodAssembler_->computeEnergy(x_) << std::endl;
//rodAssembler_->assembleGradient(x_, rhs);
rodAssembler_->assembleGradientFD(x_, rhs);
//rodAssembler_->assembleMatrix(x_, *hessianMatrix_);
rodAssembler_->assembleMatrixFD(x_, *hessianMatrix_);
//gradientFDCheck(x_, rhs, *rodAssembler_);
//hessianFDCheck(x_, *hessianMatrix_, *rodAssembler_);
rhs *= -1;
// Create trust-region obstacle on maxlevel
setTrustRegionObstacles(trustRegionRadius,
trustRegionObstacles_[grid_->maxLevel()]);
dynamic_cast<MultigridStep<MatrixType,CorrectionType>*>(mmgSolver_->iterationStep_)->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
mmgSolver_->preprocess();
dynamic_cast<MultigridStep<MatrixType,CorrectionType>*>(mmgSolver_->iterationStep_)->preprocess();
// /////////////////////////////
// Solve !
// /////////////////////////////
mmgSolver_->solve();
corr = dynamic_cast<MultigridStep<MatrixType,CorrectionType>*>(mmgSolver_->iterationStep_)->getSol();
//std::cout << "Correction: " << std::endl << corr << std::endl;
if (instrumented_) {
fprintf(fp, "Trust-region step: %d, trust-region radius: %g\n",
i, trustRegionRadius);
// ///////////////////////////////////////////////////////////////
// Compute and measure progress against the exact solution
// for each trust region step
// ///////////////////////////////////////////////////////////////
CorrectionType exactSolution = corr;
// Start from 0
double oldError = 0;
double totalConvRate = 1;
double convRate = 1;
// Write statistics of the initial solution
// Compute the energy norm
oldError = h1SemiNorm_->compute(exactSolution);
for (int j=0; j<multigridIterations_; j++) {
// read iteration from file
CorrectionType intermediateSol(grid_->size(1));
intermediateSol = 0;
char iSolFilename[100];
sprintf(iSolFilename, "tmp/mgHistory/intermediatesolution_%04d", j);
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
FILE* fpInt = fopen(iSolFilename, "rb");
if (!fpInt)
DUNE_THROW(IOError, "Couldn't open intermediate solution");
for (int k=0; k<intermediateSol.size(); k++)
for (int l=0; l<blocksize; l++)
fread(&intermediateSol[k][l], sizeof(double), 1, fpInt);
fclose(fpInt);
//std::cout << "intermediateSol\n" << intermediateSol << std::endl;
// Compute errors
intermediateSol -= exactSolution;
//std::cout << "error\n" << intermediateSol << std::endl;
// Compute the H1 norm
double error = h1SemiNorm_->compute(intermediateSol);
convRate = error / oldError;
totalConvRate *= convRate;
if (error < 1e-12)
break;
printf("Iteration: %d ", j);
printf("Errors: error %g, convergence Rate: %g, total conv rate %g\n",
error, convRate, pow(totalConvRate, 1/((double)j+1)), 0);
fprintf(fp, "%d %g %g %g\n", j+1, error, convRate, pow(totalConvRate, 1/((double)j+1)));
oldError = error;
}
}
if (this->verbosity_ == FULL)
printf("infinity norm of the correction: %g\n", corr.infinity_norm());
if (corr.infinity_norm() < tolerance_) {
if (this->verbosity_ == FULL)
std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
if (this->verbosity_ != QUIET)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
break;
}
// ////////////////////////////////////////////////////
// Check whether trust-region step can be accepted
// ////////////////////////////////////////////////////
SolutionType newIterate = x_;
for (int j=0; j<newIterate.size(); j++) {
// Add translational correction
for (int k=0; k<3; k++)
newIterate[j].r[k] += corr[j][k];
// Add rotational correction
Quaternion<double> qCorr = Quaternion<double>::exp(corr[j][3], corr[j][4], corr[j][5]);
newIterate[j].q = newIterate[j].q.mult(qCorr);
}
/** \todo Don't always recompute oldEnergy */
double oldEnergy = rodAssembler_->computeEnergy(x_);
double energy = rodAssembler_->computeEnergy(newIterate);
// compute the model decrease
// It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
// Note that rhs = -g
CorrectionType tmp(corr.size());
tmp = 0;
double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
if (this->verbosity_ == FULL) {
std::cout << "Absolute model decrease: " << modelDecrease
<< ", functional decrease: " << oldEnergy - energy << std::endl;
std::cout << "Relative model decrease: " << modelDecrease / energy
<< ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
}
assert(modelDecrease >= 0);
if (energy >= oldEnergy) {
if (this->verbosity_ == FULL)
printf("Richtung ist keine Abstiegsrichtung!\n");
#if 1
if (std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9) {
if (this->verbosity_ == FULL)
std::cout << "Suspecting rounding problems" << std::endl;
if (this->verbosity_ != QUIET)
std::cout << i+1 << " trust-region steps were taken." << std::endl;
//break;
}
#endif
// //////////////////////////////////////////////
// Check for acceptance of the step
// //////////////////////////////////////////////
if ( (oldEnergy-energy) / modelDecrease > 0.9) {
// very successful iteration
x_ = newIterate;
trustRegionRadius *= 2;
} else if ( (oldEnergy-energy) / modelDecrease > 0.01
|| std::abs(oldEnergy-energy) < 1e-12) {
// successful iteration
x_ = newIterate;
} else {
// unsuccessful iteration
trustRegionRadius /= 2;
if (this->verbosity_ == FULL)
std::cout << "Unsuccessful iteration!" << std::endl;
}
// Write current energy
if (this->verbosity_ == FULL)
std::cout << "--- Current energy: " << energy << " ---" << std::endl;
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
// /////////////////////////////////////////////////////////////////////
// Write the iterate to disk for later convergence rate measurement
// /////////////////////////////////////////////////////////////////////
if (instrumented_) {
char iRodFilename[100];
sprintf(iRodFilename, "tmp/intermediateSolution_%04d", i);
FILE* fpRod = fopen(iRodFilename, "wb");
if (!fpRod)
DUNE_THROW(SolverError, "Couldn't open file " << iRodFilename << " for writing");
for (int j=0; j<x_.size(); j++) {
for (int k=0; k<3; k++)
fwrite(&x_[j].r[k], sizeof(double), 1, fpRod);
for (int k=0; k<4; k++) // 3d hardwired here!
fwrite(&x_[j].q[k], sizeof(double), 1, fpRod);
}
fclose(fpRod);
}
// //////////////////////////////////////////////
// Close logfile
// //////////////////////////////////////////////
if (instrumented_)
fclose(fp);