diff --git a/.gitignore b/.gitignore index 5fe6fb39..5f78b452 100644 --- a/.gitignore +++ b/.gitignore @@ -27,6 +27,7 @@ linux64GccDPOpt darwinIntel64GccDPOpt # Specific files +*.log formatCode uncrustify.cfg runTests.sh diff --git a/applications/solvers/fsi/fsiFoam/fsiFoam.C b/applications/solvers/fsi/fsiFoam/fsiFoam.C index a8cbd627..d7221acf 100644 --- a/applications/solvers/fsi/fsiFoam/fsiFoam.C +++ b/applications/solvers/fsi/fsiFoam/fsiFoam.C @@ -428,16 +428,16 @@ int main( shared_ptr spaceMapping; if ( algorithm == "manifold-mapping" ) - spaceMapping = shared_ptr ( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian, initialSolutionCoarseModel ) ); + spaceMapping = shared_ptr ( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian, initialSolutionCoarseModel ) ); if ( algorithm == "output-space-mapping" ) - spaceMapping = shared_ptr ( new OutputSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); + spaceMapping = shared_ptr ( new OutputSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); if ( algorithm == "aggressive-space-mapping" ) - spaceMapping = shared_ptr ( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + spaceMapping = shared_ptr ( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); if ( algorithm == "ASM-ILS" ) - spaceMapping = shared_ptr ( new ASMILS( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + spaceMapping = shared_ptr ( new ASMILS( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); shared_ptr spaceMappingSolver( new SpaceMappingSolver( fineModel, coarseModel, spaceMapping ) ); diff --git a/src/fsi/ASMILS.C b/src/fsi/ASMILS.C index 07a1dfd1..04e32787 100644 --- a/src/fsi/ASMILS.C +++ b/src/fsi/ASMILS.C @@ -12,12 +12,13 @@ ASMILS::ASMILS( shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit ) : - SpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) + SpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) {} ASMILS::~ASMILS() @@ -49,10 +50,8 @@ void ASMILS::performPostProcessing( vector yk( m ), output( m ), R( m ), zstar( m ), zk( m ), xkprev( m ); xk = x0; xkprev = x0; - coarseResiduals.resize( m, 0 ); - coarseResiduals.setZero(); - fineResiduals.resize( m, 0 ); - fineResiduals.setZero(); + coarseResiduals.clear(); + fineResiduals.clear(); // Determine optimum of coarse model zstar if ( residualCriterium ) @@ -65,7 +64,7 @@ void ASMILS::performPostProcessing( coarseModel->optimize( y, x0, zstar ); if ( !coarseModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; + Warning << "ASMILS: surrogate model optimization process is not converged." << endl; if ( timeIndex == 0 ) xk = zstar; @@ -81,8 +80,7 @@ void ASMILS::performPostProcessing( return; } - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = xk; + fineResiduals.push_back( xk ); // Parameter extraction @@ -91,8 +89,7 @@ void ASMILS::performPostProcessing( if ( !coarseModel->allConverged() ) Warning << "Surrogate model optimization process is not converged." << endl; - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = zk; + coarseResiduals.push_back( zk ); for ( int k = 0; k < maxIter - 1; k++ ) { @@ -100,15 +97,22 @@ void ASMILS::performPostProcessing( // Determine the number of columns used to calculate the mapping matrix J - int nbCols = fineResiduals.cols() - 1; + int nbCols = fineResiduals.size() - 1; nbCols = std::max( nbCols, 0 ); - // Include information from previous time steps + // Include information from previous optimization cycles for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) - nbCols += fineResidualsList.at( i ).cols() - 1; + nbCols += fineResidualsList.at( i ).size() - 1; + + // Include information from previous time steps + + for ( unsigned i = 0; i < fineResidualsTimeList.size(); i++ ) + for ( unsigned j = 0; j < fineResidualsTimeList.at( i ).size(); j++ ) + nbCols += fineResidualsTimeList.at( i ).at( j ).size() - 1; nbCols = std::min( static_cast( xk.rows() ), nbCols ); + nbCols = std::min( nbCols, maxUsedIterations ); assert( nbCols <= xk.rows() ); @@ -125,7 +129,7 @@ void ASMILS::performPostProcessing( // Construct the V and W matrices matrix V( xk.rows(), nbCols ), W( xk.rows(), nbCols ); - int nbColsCurrentTimeStep = std::max( static_cast(fineResiduals.cols() - 1), 0 ); + int nbColsCurrentTimeStep = std::max( static_cast(fineResiduals.size() - 1), 0 ); nbColsCurrentTimeStep = std::min( nbColsCurrentTimeStep, nbCols ); // Include information from previous iterations @@ -137,29 +141,50 @@ void ASMILS::performPostProcessing( if ( colIndex >= V.cols() ) continue; - V.col( i ) = coarseResiduals.col( coarseResiduals.cols() - 2 - i ) - coarseResiduals.col( coarseResiduals.cols() - 1 ); - W.col( i ) = fineResiduals.col( fineResiduals.cols() - 2 - i ) - fineResiduals.col( fineResiduals.cols() - 1 ); + V.col( i ) = coarseResiduals.at( coarseResiduals.size() - 2 - i ) - coarseResiduals.back(); + W.col( i ) = fineResiduals.at( fineResiduals.size() - 2 - i ) - fineResiduals.back(); colIndex++; } - // Include information from previous time steps + // Include information from previous optimization cycles for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) { - assert( fineResidualsList.at( i ).cols() >= 2 ); - assert( coarseResidualsList.at( i ).cols() >= 2 ); + assert( fineResidualsList.at( i ).size() >= 2 ); + assert( coarseResidualsList.at( i ).size() >= 2 ); - for ( int j = 0; j < fineResidualsList.at( i ).cols() - 1; j++ ) + for ( unsigned j = 0; j < fineResidualsList.at( i ).size() - 1; j++ ) { if ( colIndex >= V.cols() ) continue; - V.col( colIndex ) = coarseResidualsList.at( i ).col( coarseResidualsList.at( i ).cols() - 2 - j ) - coarseResidualsList.at( i ).col( coarseResidualsList.at( i ).cols() - 1 ); - W.col( colIndex ) = fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 2 - j ) - fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 1 ); + V.col( colIndex ) = coarseResidualsList.at( i ).at( coarseResidualsList.at( i ).size() - 2 - j ) - coarseResidualsList.at( i ).back(); + W.col( colIndex ) = fineResidualsList.at( i ).at( fineResidualsList.at( i ).size() - 2 - j ) - fineResidualsList.at( i ).back(); colIndex++; } } + // Include information from previous time steps + + for ( unsigned i = 0; i < fineResidualsTimeList.size(); i++ ) + { + for ( unsigned j = 0; j < fineResidualsTimeList.at( i ).size(); j++ ) + { + assert( fineResidualsTimeList.at( i ).at( j ).size() >= 2 ); + assert( coarseResidualsTimeList.at( i ).at( j ).size() >= 2 ); + + for ( unsigned k = 0; k < fineResidualsTimeList.at( i ).at( j ).size() - 1; k++ ) + { + if ( colIndex >= V.cols() ) + continue; + + V.col( colIndex ) = coarseResidualsTimeList.at( i ).at( j ).at( coarseResidualsTimeList.at( i ).at( j ).size() - 2 - k ) - coarseResidualsTimeList.at( i ).at( j ).back(); + W.col( colIndex ) = fineResidualsTimeList.at( i ).at( j ).at( fineResidualsTimeList.at( i ).at( j ).size() - 2 - k ) - fineResidualsTimeList.at( i ).at( j ).back(); + colIndex++; + } + } + } + assert( colIndex == nbCols ); // Remove dependent columns of V and W @@ -210,28 +235,24 @@ void ASMILS::performPostProcessing( fineModel->evaluate( xk, output, R ); - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = xk; + // Check convergence criteria + if ( isConvergence( xk, xkprev, residualCriterium ) ) + { + iterationsConverged(); + break; + } + + fineResiduals.push_back( xk ); // Parameter extraction coarseModel->optimize( R, zstar, zk ); if ( !coarseModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; + Warning << "ASMILS: surrogate model optimization process is not converged." << endl; - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = zk; - - // Check convergence criteria - if ( isConvergence( xk, xkprev, residualCriterium ) ) - break; - - if ( k == maxIter - 2 ) - break; + coarseResiduals.push_back( zk ); } - - iterationsConverged(); } void ASMILS::removeColumnFromMatrix( diff --git a/src/fsi/ASMILS.H b/src/fsi/ASMILS.H index 0b8d0536..fcd6a059 100644 --- a/src/fsi/ASMILS.H +++ b/src/fsi/ASMILS.H @@ -21,6 +21,7 @@ public: shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit diff --git a/src/fsi/AggressiveSpaceMapping.C b/src/fsi/AggressiveSpaceMapping.C index 76b11465..2b0834d3 100644 --- a/src/fsi/AggressiveSpaceMapping.C +++ b/src/fsi/AggressiveSpaceMapping.C @@ -12,12 +12,13 @@ AggressiveSpaceMapping::AggressiveSpaceMapping( shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit ) : - SpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) + SpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) {} AggressiveSpaceMapping::~AggressiveSpaceMapping() @@ -38,10 +39,8 @@ void AggressiveSpaceMapping::performPostProcessing( vector yk( m ), output( m ), R( m ), zstar( m ), zk( m ), xkprev( m ); xk = x0; xkprev = x0; - coarseResiduals.resize( m, 0 ); - coarseResiduals.setZero(); - fineResiduals.resize( m, 0 ); - fineResiduals.setZero(); + coarseResiduals.clear(); + fineResiduals.clear(); // Determine optimum of coarse model zstar if ( residualCriterium ) @@ -77,11 +76,8 @@ void AggressiveSpaceMapping::performPostProcessing( if ( !coarseModel->allConverged() ) Warning << "Surrogate model optimization process is not converged." << endl; - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = zk; - - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = xk; + coarseResiduals.push_back( zk ); + fineResiduals.push_back( xk ); for ( int k = 0; k < maxIter - 1; k++ ) { @@ -93,10 +89,16 @@ void AggressiveSpaceMapping::performPostProcessing( int nbColsCurrentTimeStep = nbCols; int colIndex = 0; - // Include information from previous time steps + // Include information from previous optimization cycles for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) - nbCols += fineResidualsList.at( i ).cols() - 1; + nbCols += fineResidualsList.at( i ).size() - 1; + + // Include information from previous time steps + + for ( unsigned i = 0; i < fineResidualsTimeList.size(); i++ ) + for ( unsigned j = 0; j < fineResidualsTimeList.at( i ).size(); j++ ) + nbCols += fineResidualsTimeList.at( i ).at( j ).size() - 1; if ( nbCols > 0 ) { @@ -106,19 +108,46 @@ void AggressiveSpaceMapping::performPostProcessing( // Include information from previous time steps - for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) + for ( unsigned i = fineResidualsTimeList.size(); i-- > 0; ) { - assert( fineResidualsList.at( i ).cols() >= 2 ); - assert( coarseResidualsList.at( i ).cols() >= 2 ); + for ( unsigned j = fineResidualsTimeList.at( i ).size(); j-- > 0; ) + { + assert( fineResidualsTimeList.at( i ).at( j ).size() >= 2 ); + assert( coarseResidualsTimeList.at( i ).at( j ).size() >= 2 ); - for ( int j = 0; j < fineResidualsList.at( i ).cols() - 1; j++ ) + for ( unsigned k = 0; k < fineResidualsTimeList.at( i ).at( j ).size() - 1; k++ ) + { + colIndex++; + + vector deltax, deltaz; + + deltax = fineResidualsTimeList.at( i ).at( j ).at( k + 1 ) - fineResidualsTimeList.at( i ).at( j ).at( k ); + deltaz = coarseResidualsTimeList.at( i ).at( j ).at( k + 1 ) - coarseResidualsTimeList.at( i ).at( j ).at( k ); + + if ( deltax.norm() >= singularityLimit ) + { + // Sherman–Morrison formula + J += (deltax - J * deltaz) / (deltax.transpose() * J * deltaz) * (deltax.transpose() * J); + } + } + } + } + + // Include information from previous optimization cycles + + for ( unsigned i = fineResidualsList.size(); i-- > 0; ) + { + assert( fineResidualsList.at( i ).size() >= 2 ); + assert( coarseResidualsList.at( i ).size() >= 2 ); + + for ( unsigned j = 0; j < fineResidualsList.at( i ).size() - 1; j++ ) { colIndex++; vector deltax, deltaz; - deltax = fineResidualsList.at( i ).col( j + 1 ) - fineResidualsList.at( i ).col( j ); - deltaz = coarseResidualsList.at( i ).col( j + 1 ) - coarseResidualsList.at( i ).col( j ); + deltax = fineResidualsList.at( i ).at( j + 1 ) - fineResidualsList.at( i ).at( j ); + deltaz = coarseResidualsList.at( i ).at( j + 1 ) - coarseResidualsList.at( i ).at( j ); if ( deltax.norm() >= singularityLimit ) { @@ -136,8 +165,8 @@ void AggressiveSpaceMapping::performPostProcessing( vector deltax, deltaz; - deltax = fineResiduals.col( i + 1 ) - fineResiduals.col( i ); - deltaz = coarseResiduals.col( i + 1 ) - coarseResiduals.col( i ); + deltax = fineResiduals.at( i + 1 ) - fineResiduals.at( i ); + deltaz = coarseResiduals.at( i + 1 ) - coarseResiduals.at( i ); if ( deltax.norm() >= singularityLimit ) { @@ -159,7 +188,10 @@ void AggressiveSpaceMapping::performPostProcessing( // Check convergence criteria if ( isConvergence( xk, xkprev, residualCriterium ) ) + { + iterationsConverged(); break; + } // Parameter extraction @@ -168,12 +200,7 @@ void AggressiveSpaceMapping::performPostProcessing( if ( !coarseModel->allConverged() ) Warning << "Surrogate model optimization process is not converged." << endl; - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = zk; - - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = xk; + coarseResiduals.push_back( zk ); + fineResiduals.push_back( xk ); } - - iterationsConverged(); } diff --git a/src/fsi/AggressiveSpaceMapping.H b/src/fsi/AggressiveSpaceMapping.H index 240fce99..ab3c4cd5 100644 --- a/src/fsi/AggressiveSpaceMapping.H +++ b/src/fsi/AggressiveSpaceMapping.H @@ -21,6 +21,7 @@ public: shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit diff --git a/src/fsi/AitkenPostProcessing.C b/src/fsi/AitkenPostProcessing.C index 2f3768a2..26704a66 100644 --- a/src/fsi/AitkenPostProcessing.C +++ b/src/fsi/AitkenPostProcessing.C @@ -119,11 +119,12 @@ void AitkenPostProcessing::performPostProcessing( // Check convergence criteria if ( isConvergence( xk, xkprev, residualCriterium ) ) + { + bool keepIterations = residualCriterium; + iterationsConverged( keepIterations ); break; + } assert( fsi->iter <= maxIter ); } - - bool keepIterations = residualCriterium; - iterationsConverged( keepIterations ); } diff --git a/src/fsi/AndersonPostProcessing.C b/src/fsi/AndersonPostProcessing.C index cc9c0cbd..8571534b 100644 --- a/src/fsi/AndersonPostProcessing.C +++ b/src/fsi/AndersonPostProcessing.C @@ -368,7 +368,15 @@ namespace fsi // Check convergence criteria if ( isConvergence( output, output + y - R, residualCriterium ) ) + { + bool keepIterations = residualCriterium || solsList.size() == 0; + iterationsConverged( keepIterations ); + + if ( updateJacobian && J.cols() > 0 && timeIndex >= reuseInformationStartingFromTimeIndex ) + Jprev = J; + break; + } determineScalingFactors( output ); @@ -383,12 +391,6 @@ namespace fsi assert( sols.at( 0 ).rows() == residuals.at( 0 ).rows() ); assert( fsi->iter <= maxIter ); } - - bool keepIterations = residualCriterium || solsList.size() == 0; - iterationsConverged( keepIterations ); - - if ( updateJacobian && J.cols() > 0 && timeIndex >= reuseInformationStartingFromTimeIndex ) - Jprev = J; } } diff --git a/src/fsi/BroydenPostProcessing.C b/src/fsi/BroydenPostProcessing.C index 63bd7bd1..7810dea5 100644 --- a/src/fsi/BroydenPostProcessing.C +++ b/src/fsi/BroydenPostProcessing.C @@ -249,13 +249,14 @@ void BroydenPostProcessing::performPostProcessing( // Check convergence criteria if ( isConvergence( output, output + y - R, residualCriterium ) ) + { + bool keepIterations = residualCriterium; + iterationsConverged( keepIterations ); break; + } assert( sols.size() == residuals.size() ); assert( sols.at( 0 ).rows() == residuals.at( 0 ).rows() ); assert( fsi->iter <= maxIter ); } - - bool keepIterations = residualCriterium; - iterationsConverged( keepIterations ); } diff --git a/src/fsi/ManifoldMapping.C b/src/fsi/ManifoldMapping.C index 0b8bead1..fd46c3d3 100644 --- a/src/fsi/ManifoldMapping.C +++ b/src/fsi/ManifoldMapping.C @@ -8,17 +8,34 @@ using namespace fsi; +ManifoldMapping::ManifoldMapping( + shared_ptr fineModel, + shared_ptr coarseModel, + int maxIter + ) + : + SpaceMapping( fineModel, coarseModel, maxIter, 500, 0, 0, 1.0e-15 ), + updateJacobian( false ), + initialSolutionCoarseModel( false ), + scaling( false ), + iter( 0 ), + scalingFactors( Eigen::VectorXd::Ones( 2 ) ), + sizeVar0( 0 ), + sizeVar1( 0 ) +{} + ManifoldMapping::ManifoldMapping( shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit, bool updateJacobian ) : - SpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), + SpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), updateJacobian( updateJacobian ), initialSolutionCoarseModel( false ), scaling( false ), @@ -32,6 +49,7 @@ ManifoldMapping::ManifoldMapping( shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit, @@ -39,7 +57,7 @@ ManifoldMapping::ManifoldMapping( bool initialSolutionCoarseModel ) : - SpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), + SpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), updateJacobian( updateJacobian ), initialSolutionCoarseModel( initialSolutionCoarseModel ), scaling( false ), @@ -113,10 +131,8 @@ void ManifoldMapping::performPostProcessing( vector xkprev = x0; output.setZero(); R.setZero(); - coarseResiduals.resize( m, 1 ); - coarseResiduals.setZero(); - fineResiduals.resize( m, 1 ); - fineResiduals.setZero(); + coarseResiduals.clear(); + fineResiduals.clear(); iter = 0; matrix Tk; @@ -144,7 +160,7 @@ void ManifoldMapping::performPostProcessing( // Coarse model evaluation coarseModel->evaluate( xk, output, R ); - coarseResiduals.col( 0 ) = R; + coarseResiduals.push_back( R ); assert( xk.rows() == n ); assert( output.rows() == m ); assert( R.rows() == m ); @@ -152,7 +168,7 @@ void ManifoldMapping::performPostProcessing( // Fine model evaluation fineModel->evaluate( xk, output, R ); - fineResiduals.col( 0 ) = R; + fineResiduals.push_back( R ); assert( output.rows() == m ); assert( R.rows() == m ); @@ -179,17 +195,25 @@ void ManifoldMapping::performPostProcessing( int nbCols = std::min( k, n ); int nbColsCurrentTimeStep = nbCols; - // Include information from previous time steps + // Include information from previous optimization cycles for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) - nbCols += fineResidualsList.at( i ).cols() - 1; + nbCols += fineResidualsList.at( i ).size() - 1; + + // Include information from previous time steps + + for ( unsigned i = 0; i < fineResidualsTimeList.size(); i++ ) + for ( unsigned j = 0; j < fineResidualsTimeList.at( i ).size(); j++ ) + nbCols += fineResidualsTimeList.at( i ).at( j ).size() - 1; nbCols = std::min( nbCols, n ); + nbCols = std::min( nbCols, maxUsedIterations ); + nbColsCurrentTimeStep = std::min( nbCols, nbColsCurrentTimeStep ); // Update the design specification yk - vector alpha = fineResiduals.col( k ) - y; - yk = coarseResiduals.col( k ); + vector alpha = fineResiduals.at( k ) - y; + yk = coarseResiduals.at( k ); // Apply scaling if ( scaling ) @@ -216,30 +240,52 @@ void ManifoldMapping::performPostProcessing( int colIndex = 0; + // Include information from previous iterations + for ( int i = 0; i < nbColsCurrentTimeStep; i++ ) { - DeltaF.col( i ) = fineResiduals.col( k ) - fineResiduals.col( k - 1 - i ); - DeltaC.col( i ) = coarseResiduals.col( k ) - coarseResiduals.col( k - 1 - i ); + DeltaF.col( i ) = fineResiduals.back() - fineResiduals.at( k - 1 - i ); + DeltaC.col( i ) = coarseResiduals.back() - coarseResiduals.at( k - 1 - i ); colIndex++; } - // Include information from previous time steps + // Include information from previous optimization cycles for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) { - assert( fineResidualsList.at( i ).cols() >= 2 ); + assert( fineResidualsList.at( i ).size() >= 2 ); - for ( unsigned j = 0; j < fineResidualsList.at( i ).cols() - 1; j++ ) + for ( unsigned j = 0; j < fineResidualsList.at( i ).size() - 1; j++ ) { if ( colIndex >= DeltaF.cols() ) continue; - DeltaF.col( colIndex ) = fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 1 ) - fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 2 - j ); - DeltaC.col( colIndex ) = coarseResidualsList.at( i ).col( coarseResidualsList.at( i ).cols() - 1 ) - coarseResidualsList.at( i ).col( coarseResidualsList.at( i ).cols() - 2 - j ); + DeltaF.col( colIndex ) = fineResidualsList.at( i ).back() - fineResidualsList.at( i ).at( fineResidualsList.at( i ).size() - 2 - j ); + DeltaC.col( colIndex ) = coarseResidualsList.at( i ).back() - coarseResidualsList.at( i ).at( coarseResidualsList.at( i ).size() - 2 - j ); colIndex++; } } + // Include information from previous time steps + + for ( unsigned i = 0; i < fineResidualsTimeList.size(); i++ ) + { + for ( unsigned j = 0; j < fineResidualsTimeList.at( i ).size(); j++ ) + { + assert( fineResidualsTimeList.at( i ).at( j ).size() >= 2 ); + + for ( unsigned k = 0; k < fineResidualsTimeList.at( i ).at( j ).size() - 1; k++ ) + { + if ( colIndex >= DeltaF.cols() ) + continue; + + DeltaF.col( colIndex ) = fineResidualsTimeList.at( i ).at( j ).back() - fineResidualsTimeList.at( i ).at( j ).at( fineResidualsTimeList.at( i ).at( j ).size() - 2 - k ); + DeltaC.col( colIndex ) = coarseResidualsTimeList.at( i ).at( j ).back() - coarseResidualsTimeList.at( i ).at( j ).at( coarseResidualsTimeList.at( i ).at( j ).size() - 2 - k ); + colIndex++; + } + } + } + assert( colIndex == nbCols ); // Apply scaling @@ -369,8 +415,7 @@ void ManifoldMapping::performPostProcessing( assert( xk.rows() == n ); assert( output.rows() == m ); assert( R.rows() == m ); - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = R; + coarseResiduals.push_back( R ); // Fine model evaluation @@ -378,8 +423,7 @@ void ManifoldMapping::performPostProcessing( assert( xk.rows() == n ); assert( output.rows() == m ); assert( R.rows() == m ); - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = R; + fineResiduals.push_back( R ); determineScalingFactors( output ); @@ -393,11 +437,11 @@ void ManifoldMapping::performPostProcessing( if ( updateJacobian && Tk.cols() > 0 && timeIndex >= reuseInformationStartingFromTimeIndex ) Tkprev = Tk; + iterationsConverged(); + break; } } - - iterationsConverged(); } void ManifoldMapping::removeColumnFromMatrix( diff --git a/src/fsi/ManifoldMapping.H b/src/fsi/ManifoldMapping.H index bae38010..21fa5371 100644 --- a/src/fsi/ManifoldMapping.H +++ b/src/fsi/ManifoldMapping.H @@ -18,10 +18,17 @@ public: using SpaceMapping::performPostProcessing; + ManifoldMapping( + shared_ptr fineModel, + shared_ptr coarseModel, + int maxIter + ); + ManifoldMapping( shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit, @@ -32,6 +39,7 @@ public: shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit, diff --git a/src/fsi/OutputSpaceMapping.C b/src/fsi/OutputSpaceMapping.C index 4cee2566..786b11ad 100644 --- a/src/fsi/OutputSpaceMapping.C +++ b/src/fsi/OutputSpaceMapping.C @@ -12,13 +12,14 @@ OutputSpaceMapping::OutputSpaceMapping( shared_ptr fineModel, shared_ptr surrogateModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit, int order ) : - SpaceMapping( fineModel, surrogateModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), + SpaceMapping( fineModel, surrogateModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), surrogateModel( surrogateModel ), order( order ) { @@ -29,23 +30,37 @@ OutputSpaceMapping::OutputSpaceMapping( OutputSpaceMapping::~OutputSpaceMapping() {} -void OutputSpaceMapping::iterationsConverged() +void OutputSpaceMapping::finalizeTimeStep() { - SpaceMapping::iterationsConverged(); + SpaceMapping::finalizeTimeStep(); // Save input/output information for next time step - if ( nbReuse > 0 && sols.size() >= 2 && timeIndex >= reuseInformationStartingFromTimeIndex ) - solsList.push_front( sols ); + if ( nbReuse > 0 && solsList.size() >= 1 && timeIndex >= reuseInformationStartingFromTimeIndex ) + solsTimeList.push_front( solsList ); - // Remove the last item from the sols list and solutions - // list in order to ensure that at maximum nbReuse time steps + // Remove the last items from the residual list and solutions list + // in order to ensure that at maximum nbReuse time steps // are included. - while ( static_cast( solsList.size() ) > nbReuse ) - solsList.pop_back(); + while ( static_cast( solsTimeList.size() ) > nbReuse ) + solsTimeList.pop_back(); + + solsList.clear(); + + assert( solsTimeList.size() == coarseResidualsTimeList.size() ); + assert( solsTimeList.size() == fineResidualsTimeList.size() ); + assert( static_cast( solsTimeList.size() ) <= nbReuse ); +} + +void OutputSpaceMapping::iterationsConverged() +{ + SpaceMapping::iterationsConverged(); + + // Save input/output information for next solve + if ( sols.size() >= 2 ) + solsList.push_front( sols ); sols.clear(); - assert( static_cast( solsList.size() ) <= nbReuse ); assert( solsList.size() == coarseResidualsList.size() ); assert( solsList.size() == fineResidualsList.size() ); } @@ -68,10 +83,8 @@ void OutputSpaceMapping::performPostProcessing( vector xkprev = x0; output.setZero(); R.setZero(); - coarseResiduals.resize( m, 1 ); - coarseResiduals.setZero(); - fineResiduals.resize( m, 1 ); - fineResiduals.setZero(); + coarseResiduals.clear(); + fineResiduals.clear(); if ( timeIndex == 0 ) { @@ -86,7 +99,7 @@ void OutputSpaceMapping::performPostProcessing( surrogateModel->optimize( y, x0, xk ); if ( !surrogateModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; + Warning << "Output space mapping: surrogate model optimization process is not converged." << endl; } assert( xk.rows() == n ); @@ -99,7 +112,7 @@ void OutputSpaceMapping::performPostProcessing( // Coarse model evaluation surrogateModel->evaluate( xk, output, R ); - coarseResiduals.col( 0 ) = R; + coarseResiduals.push_back( R ); assert( xk.rows() == n ); assert( output.rows() == m ); assert( R.rows() == m ); @@ -107,7 +120,7 @@ void OutputSpaceMapping::performPostProcessing( // Fine model evaluation fineModel->evaluate( xk, output, R ); - fineResiduals.col( 0 ) = R; + fineResiduals.push_back( R ); assert( output.rows() == m ); assert( R.rows() == m ); @@ -130,18 +143,23 @@ void OutputSpaceMapping::performPostProcessing( int nbCols = k; int nbColsCurrentTimeStep = nbCols; - // Include information from previous time steps + // Include information from previous optimization cycles + + for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) + nbCols += fineResidualsList.at( i ).size() - 1; - for ( unsigned i = 0; i < solsList.size(); i++ ) - nbCols += solsList.at( i ).size() - 1; + // Include information from previous time steps + for ( unsigned i = 0; i < fineResidualsTimeList.size(); i++ ) + for ( unsigned j = 0; j < fineResidualsTimeList.at( i ).size(); j++ ) + nbCols += fineResidualsTimeList.at( i ).at( j ).size() - 1; // Update the design specification yk - yk = coarseResiduals.col( k ) - (fineResiduals.col( k ) - y); + yk = coarseResiduals.at( k ) - (fineResiduals.at( k ) - y); if ( nbCols > 0 && order == 1 ) { - assert( fineResiduals.cols() == coarseResiduals.cols() ); - assert( static_cast( sols.size() ) == fineResiduals.cols() ); + assert( fineResiduals.size() == coarseResiduals.size() ); + assert( sols.size() == fineResiduals.size() ); assert( solsList.size() == fineResidualsList.size() ); assert( solsList.size() == coarseResidualsList.size() ); @@ -150,20 +168,45 @@ void OutputSpaceMapping::performPostProcessing( fsi::vector d, dprev, deltad, deltax; int colIndex = 0; - Info << "Output space mapping with "; + Info << "OSM(1) with "; Info << nbCols; Info << " cols for the Jacobian" << endl; // Include information from previous time steps + for ( unsigned i = solsTimeList.size(); i-- > 0; ) + { + for ( unsigned j = solsTimeList.at( i ).size(); j-- > 0; ) + { + for ( unsigned k = 0; k < solsTimeList.at( i ).at( j ).size() - 1; k++ ) + { + colIndex++; + + d = fineResidualsTimeList.at( i ).at( j ).at( k + 1 ) - coarseResidualsTimeList.at( i ).at( j ).at( k + 1 ); + dprev = fineResidualsTimeList.at( i ).at( j ).at( k ) - coarseResidualsTimeList.at( i ).at( j ).at( k ); + deltad = d - dprev; + deltax = solsTimeList.at( i ).at( j ).at( k + 1 ) - solsTimeList.at( i ).at( j ).at( k ); + + // Broyden update for the Jacobian matrix + + if ( deltax.norm() < singularityLimit ) + continue; + + J += (deltad - J * deltax) / deltax.squaredNorm() * deltax.transpose(); + } + } + } + + // Include information from previous optimization cycles + for ( unsigned i = solsList.size(); i-- > 0; ) { for ( unsigned j = 0; j < solsList.at( i ).size() - 1; j++ ) { colIndex++; - d = fineResidualsList.at( i ).col( j + 1 ) - coarseResidualsList.at( i ).col( j + 1 ); - dprev = fineResidualsList.at( i ).col( j ) - coarseResidualsList.at( i ).col( j ); + d = fineResidualsList.at( i ).at( j + 1 ) - coarseResidualsList.at( i ).at( j + 1 ); + dprev = fineResidualsList.at( i ).at( j ) - coarseResidualsList.at( i ).at( j ); deltad = d - dprev; deltax = solsList.at( i ).at( j + 1 ) - solsList.at( i ).at( j ); @@ -180,8 +223,8 @@ void OutputSpaceMapping::performPostProcessing( { colIndex++; - d = fineResiduals.col( i + 1 ) - coarseResiduals.col( i + 1 ); - dprev = fineResiduals.col( i ) - coarseResiduals.col( i ); + d = fineResiduals.at( i + 1 ) - coarseResiduals.at( i + 1 ); + dprev = fineResiduals.at( i ) - coarseResiduals.at( i ); deltad = d - dprev; deltax = sols.at( i + 1 ) - sols.at( i ); @@ -202,6 +245,12 @@ void OutputSpaceMapping::performPostProcessing( if ( nbCols > 0 && order == 2 ) { nbCols = std::min( nbCols, n ); + nbCols = std::min( nbCols, maxUsedIterations ); + nbColsCurrentTimeStep = std::min( nbCols, nbColsCurrentTimeStep ); + + Info << "OSM(2) with "; + Info << nbCols; + Info << " cols for the Jacobian" << endl; matrix DeltaF( m, nbCols ), DeltaX( m, nbCols ); @@ -214,13 +263,13 @@ void OutputSpaceMapping::performPostProcessing( if ( colIndex >= DeltaF.cols() ) continue; - DeltaF.col( colIndex ) = fineResiduals.col( k ) - coarseResiduals.col( k ); - DeltaF.col( colIndex ) -= fineResiduals.col( k - 1 - i ) - coarseResiduals.col( k - 1 - i ); - DeltaX.col( colIndex ) = sols.at( k ) - sols.at( k - 1 - i ); + DeltaF.col( colIndex ) = fineResiduals.back() - coarseResiduals.back(); + DeltaF.col( colIndex ) -= fineResiduals.at( k - 1 - i ) - coarseResiduals.at( k - 1 - i ); + DeltaX.col( colIndex ) = sols.back() - sols.at( k - 1 - i ); colIndex++; } - // Include information from previous time steps + // Include information from previous optimization cycles for ( unsigned i = 0; i < solsList.size(); i++ ) { @@ -229,16 +278,44 @@ void OutputSpaceMapping::performPostProcessing( if ( colIndex >= DeltaF.cols() ) continue; - DeltaF.col( colIndex ) = fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 1 ) - coarseResidualsList.at( i ).col( coarseResidualsList.at( i ).cols() - 1 ); + DeltaF.col( colIndex ) = fineResidualsList.at( i ).back() - coarseResidualsList.at( i ).back(); - DeltaF.col( colIndex ) -= fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 2 - j ) - coarseResidualsList.at( i ).col( coarseResidualsList.at( i ).cols() - 2 - j ); + DeltaF.col( colIndex ) -= fineResidualsList.at( i ).at( fineResidualsList.at( i ).size() - 2 - j ); + DeltaF.col( colIndex ) -= coarseResidualsList.at( i ).at( coarseResidualsList.at( i ).size() - 2 - j ); - DeltaX.col( colIndex ) = solsList.at( i ).at( solsList.at( i ).size() - 1 ) - solsList.at( i ).at( solsList.at( i ).size() - 2 - j ); + DeltaX.col( colIndex ) = solsList.at( i ).back() - solsList.at( i ).at( solsList.at( i ).size() - 2 - j ); colIndex++; } } + // Include information from previous time steps + for ( unsigned i = 0; i < solsTimeList.size(); i++ ) + { + for ( unsigned j = 0; j < solsTimeList.at( i ).size(); j++ ) + { + for ( unsigned k = 0; k < solsTimeList.at( i ).at( j ).size() - 1; k++ ) + { + assert( fineResidualsTimeList.at( i ).at( j ).size() == solsTimeList.at( i ).at( j ).size() ); + assert( coarseResidualsTimeList.at( i ).at( j ).size() == solsTimeList.at( i ).at( j ).size() ); + + if ( colIndex >= DeltaF.cols() ) + continue; + + DeltaF.col( colIndex ) = fineResidualsTimeList.at( i ).at( j ).back(); + DeltaF.col( colIndex ) -= coarseResidualsTimeList.at( i ).at( j ).back(); + + DeltaF.col( colIndex ) -= fineResidualsTimeList.at( i ).at( j ).at( fineResidualsTimeList.at( i ).at( j ).size() - 2 - k ); + DeltaF.col( colIndex ) -= coarseResidualsTimeList.at( i ).at( j ).at( coarseResidualsTimeList.at( i ).at( j ).size() - 2 - k ); + + DeltaX.col( colIndex ) = solsTimeList.at( i ).at( j ).back(); + DeltaX.col( colIndex ) -= solsTimeList.at( i ).at( j ).at( solsTimeList.at( i ).at( j ).size() - 2 - k ); + + colIndex++; + } + } + } + assert( colIndex == nbCols ); // Remove dependent columns of DeltaC and DeltaF @@ -296,7 +373,7 @@ void OutputSpaceMapping::performPostProcessing( assert( xk.rows() == n ); if ( !surrogateModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; + Warning << "Output space mapping: surrogate model optimization process is not converged." << endl; xk = output; @@ -308,8 +385,7 @@ void OutputSpaceMapping::performPostProcessing( assert( xk.rows() == n ); assert( output.rows() == m ); assert( R.rows() == m ); - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = R; + coarseResiduals.push_back( R ); // Fine model evaluation @@ -317,8 +393,7 @@ void OutputSpaceMapping::performPostProcessing( assert( xk.rows() == n ); assert( output.rows() == m ); assert( R.rows() == m ); - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = R; + fineResiduals.push_back( R ); sols.push_back( xk ); @@ -326,11 +401,10 @@ void OutputSpaceMapping::performPostProcessing( if ( isConvergence( xk, xkprev, residualCriterium ) ) { assert( fineModel->allConverged() ); + iterationsConverged(); break; } } - - iterationsConverged(); } void OutputSpaceMapping::removeColumnFromMatrix( diff --git a/src/fsi/OutputSpaceMapping.H b/src/fsi/OutputSpaceMapping.H index 33f7e6f2..89d1ce5d 100644 --- a/src/fsi/OutputSpaceMapping.H +++ b/src/fsi/OutputSpaceMapping.H @@ -26,6 +26,7 @@ public: shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit, @@ -34,6 +35,8 @@ public: virtual ~OutputSpaceMapping(); + virtual void finalizeTimeStep(); + virtual void iterationsConverged(); virtual void performPostProcessing( @@ -51,6 +54,7 @@ public: shared_ptr surrogateModel; deque sols; deque > solsList; + deque > > solsTimeList; int order; }; diff --git a/src/fsi/PostProcessing.H b/src/fsi/PostProcessing.H index 3f239e1a..2775cbec 100644 --- a/src/fsi/PostProcessing.H +++ b/src/fsi/PostProcessing.H @@ -9,7 +9,6 @@ #include "MultiLevelFsiSolver.H" -using std::list; using std::deque; namespace fsi diff --git a/src/fsi/ResponseParameterMapping.C b/src/fsi/ResponseParameterMapping.C deleted file mode 100644 index f44a7cdb..00000000 --- a/src/fsi/ResponseParameterMapping.C +++ /dev/null @@ -1,421 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#include "ResponseParameterMapping.H" - -using namespace fsi; - -ResponseParameterMapping::ResponseParameterMapping( - shared_ptr fineModel, - shared_ptr coarseModel, - int maxIter, - int nbReuse, - int reuseInformationStartingFromTimeIndex, - double singularityLimit - ) - : - SpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ), - iter( 0 ) -{} - -ResponseParameterMapping::~ResponseParameterMapping() -{} - -void ResponseParameterMapping::iterationsConverged() -{ - SpaceMapping::iterationsConverged(); - - assert( xktildeList.size() == xkList.size() ); - - // Save input/output information for next time step - if ( nbReuse > 0 && xktildeList.size() >= 2 && timeIndex >= reuseInformationStartingFromTimeIndex ) - { - xktildePreviousTimeStepsList.push_front( xktildeList ); - xkPreviousTimeStepsList.push_front( xkList ); - } - - // Remove the last item from the sols list and solutions - // list in order to ensure that at maximum nbReuse time steps - // are included. - while ( static_cast( xktildePreviousTimeStepsList.size() ) > nbReuse ) - { - xktildePreviousTimeStepsList.pop_back(); - xkPreviousTimeStepsList.pop_back(); - } - - xktildeList.clear(); - xkList.clear(); - - assert( static_cast( xktildePreviousTimeStepsList.size() ) <= nbReuse ); - assert( xkPreviousTimeStepsList.size() == coarseResidualsList.size() ); - assert( xktildePreviousTimeStepsList.size() == fineResidualsList.size() ); - assert( xktildePreviousTimeStepsList.size() == xkPreviousTimeStepsList.size() ); -} - -void ResponseParameterMapping::performPostProcessing( - const vector & y, - const vector & x0, - vector & xk, - bool residualCriterium - ) -{ - assert( x0.rows() == xk.rows() ); - - // Initialize variables - - int m = y.rows(); - int n = x0.rows(); - matrix B( x0.rows(), x0.rows() ); - B.setZero(); - vector yk( m ), output( m ), R( m ); - vector xktilde = x0; - xk = x0; - vector xkprev = x0; - output.setZero(); - R.setZero(); - coarseResiduals.resize( m, 1 ); - coarseResiduals.setZero(); - fineResiduals.resize( m, 1 ); - fineResiduals.setZero(); - iter = 0; - xktildeList.clear(); - xkList.clear(); - mappedCoarseList.clear(); - - // Determine optimum of coarse model xstar - if ( residualCriterium ) - { - assert( y.norm() < 1.0e-14 ); - coarseModel->optimize( x0, xk ); - } - - if ( !residualCriterium ) - coarseModel->optimize( y, x0, xk ); - - if ( !coarseModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; - - // Initialize coarse model and fine model responses - - // Fine model evaluation - - fineModel->evaluate( xk, output, R ); - fineResiduals.col( 0 ) = R; - - iter++; - - // Parameter extraction - Info << nl << "Response and parameter mapping: parameter extraction process" << endl; - coarseModel->optimize( R, xk, xktilde ); - xktildeList.push_back( xktilde ); - xkList.push_back( xk ); - - // Coarse model evaluation - - coarseModel->evaluate( xktilde, output, R ); - coarseResiduals.col( 0 ) = R; - - // Check convergence criteria - if ( isConvergence( xk, xkprev, residualCriterium ) ) - { - iterationsConverged(); - return; - } - - for ( int k = 0; k < maxIter - 1; k++ ) - { - xkprev = xk; - - // Determine the number of columns used to calculate the mapping matrix - - int nbCols = std::min( k, n ); - int nbColsCurrentTimeStep = nbCols; - - // Include information from previous time steps - - for ( std::deque::iterator it = fineResidualsList.begin(); it != fineResidualsList.end(); ++it ) - nbCols += it->cols() - 1; - - nbCols = std::min( nbCols, n ); - - if ( nbCols == 0 ) - // Update the design specification yk - yk = coarseResiduals.col( k ) - (fineResiduals.col( k ) - y); - - if ( nbCols >= 1 ) - { - // Update the mapping matrix B - - matrix DeltaXTilde( m, nbCols ), DeltaX( m, nbCols ); - - int colIndex = 0; - - for ( int i = 0; i < nbColsCurrentTimeStep; i++ ) - { - DeltaXTilde.col( i ) = xktildeList.at( k - 1 - i ) - xktildeList.at( k ); - DeltaX.col( i ) = xkList.at( k - 1 - i ) - xkList.at( k ); - colIndex++; - } - - // Include information from previous time steps for the mapping matrix B - - assert( xkPreviousTimeStepsList.size() == xktildePreviousTimeStepsList.size() ); - assert( static_cast ( fineResidualsList.size() ) == static_cast ( xkPreviousTimeStepsList.size() ) ); - - for ( unsigned i = 0; i < xkPreviousTimeStepsList.size(); i++ ) - { - assert( xkPreviousTimeStepsList.at( i ).size() == xktildePreviousTimeStepsList.at( i ).size() ); - assert( static_cast( xkPreviousTimeStepsList.at( i ).size() ) == fineResidualsList.at( i ).cols() ); - - for ( unsigned j = 0; j < xkPreviousTimeStepsList.at( i ).size() - 1; j++ ) - { - if ( colIndex >= DeltaXTilde.cols() ) - continue; - - DeltaXTilde.col( colIndex ) = xktildePreviousTimeStepsList.at( i ).at( xktildePreviousTimeStepsList.at( i ).size() - 2 - j ) - xktildePreviousTimeStepsList.at( i ).at( xktildePreviousTimeStepsList.at( i ).size() - 1 ); - DeltaX.col( colIndex ) = xkPreviousTimeStepsList.at( i ).at( xkPreviousTimeStepsList.at( i ).size() - 2 - j ) - xkPreviousTimeStepsList.at( i ).at( xkPreviousTimeStepsList.at( i ).size() - 1 ); - colIndex++; - } - } - - assert( colIndex == nbCols ); - - // Remove dependent columns of DeltaX and DeltaXTilde - - int nbRemoveCols = 1; - - while ( nbRemoveCols > 0 ) - { - nbRemoveCols = 0; - - if ( DeltaX.cols() == 1 ) - break; - - // Calculate singular value decomposition with Eigen - - Eigen::JacobiSVD svd( DeltaX, Eigen::ComputeThinU | Eigen::ComputeThinV ); - - vector singularValues = svd.singularValues(); - - for ( int i = 0; i < singularValues.rows(); i++ ) - { - if ( std::abs( singularValues( i ) ) <= singularityLimit ) - { - // Remove the column from DeltaC and DeltaF - removeColumnFromMatrix( DeltaX, i - nbRemoveCols ); - removeColumnFromMatrix( DeltaXTilde, i - nbRemoveCols ); - - nbRemoveCols++; - } - } - - if ( nbRemoveCols ) - Info << "Response and parameter mapping: remove " << nbRemoveCols << " columns from the DeltaX and DeltaXTilde matrices" << endl; - } - - assert( DeltaX.cols() == DeltaXTilde.cols() ); - assert( DeltaX.cols() >= 1 ); - - // Calculate singular value decomposition with Eigen - - Eigen::JacobiSVD svd_X( DeltaX, Eigen::ComputeThinU | Eigen::ComputeThinV ); - Eigen::JacobiSVD svd_XTilde( DeltaXTilde, Eigen::ComputeThinU | Eigen::ComputeThinV ); - - matrix pseudoSigma_X = svd_X.singularValues().asDiagonal(); - - for ( int i = 0; i < pseudoSigma_X.cols(); i++ ) - pseudoSigma_X( i, i ) = 1.0 / pseudoSigma_X( i, i ); - - matrix pseudoDeltaX = svd_X.matrixV() * pseudoSigma_X * svd_X.matrixU().transpose(); - matrix I = Eigen::MatrixXd::Identity( m, m ); - - // Calculate the new mapping matrix B - - B = DeltaXTilde * pseudoDeltaX + ( I - svd_XTilde.matrixU() * svd_XTilde.matrixU().transpose() ) * ( I - svd_X.matrixU() * svd_X.matrixU().transpose() ); - - // Evaluate the coarse model for each xk to fill the mapped coarse response list mappedCoarseList - mappedCoarseList.clear(); - vector xkTmp, outputTmp( m ), Rtmp( m ); - - for ( unsigned i = 0; i < xkList.size(); i++ ) - { - xkTmp = xkList.at( i ); - coarseModel->evaluate( xktilde + B * (xkTmp - xk), outputTmp, Rtmp ); - mappedCoarseList.push_back( Rtmp ); - } - - assert( mappedCoarseList.size() == xkList.size() ); - - // Evaluate the coarse model for each xk for the previous time steps - - mappedCoarsePreviousTimeStepsList.clear(); - - for ( unsigned i = 0; i < xkPreviousTimeStepsList.size(); i++ ) - { - std::deque mappedCoarsePreviousTimeStepList; - - for ( unsigned j = 0; j < xkPreviousTimeStepsList.at( i ).size(); j++ ) - { - xkTmp = xkPreviousTimeStepsList.at( i ).at( j ); - coarseModel->evaluate( xktilde + B * (xkTmp - xk), outputTmp, Rtmp ); - mappedCoarsePreviousTimeStepList.push_back( Rtmp ); - } - - mappedCoarsePreviousTimeStepsList.push_back( mappedCoarsePreviousTimeStepList ); - - assert( mappedCoarsePreviousTimeStepList.size() == xkPreviousTimeStepsList.at( i ).size() ); - } - - assert( mappedCoarsePreviousTimeStepsList.size() == xkPreviousTimeStepsList.size() ); - - // Update the mapping matrix T - - matrix DeltaF( m, nbCols ), DeltaC( m, nbCols ); - - colIndex = 0; - - for ( int i = 0; i < nbColsCurrentTimeStep; i++ ) - { - DeltaF.col( i ) = fineResiduals.col( k - 1 - i ) - fineResiduals.col( k ); - DeltaC.col( i ) = mappedCoarseList.at( k - 1 - i ) - mappedCoarseList.at( k ); - colIndex++; - } - - // Include information from previous time steps for the mapping matrix T - - for ( unsigned i = 0; i < fineResidualsList.size(); i++ ) - { - for ( unsigned j = 0; j < fineResidualsList.at( i ).cols() - 1; j++ ) - { - if ( colIndex >= DeltaF.cols() ) - continue; - - DeltaF.col( colIndex ) = fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 2 - j ) - fineResidualsList.at( i ).col( fineResidualsList.at( i ).cols() - 1 ); - DeltaC.col( colIndex ) = mappedCoarsePreviousTimeStepsList.at( i ).at( mappedCoarsePreviousTimeStepsList.at( i ).size() - 2 - j ) - mappedCoarsePreviousTimeStepsList.at( i ).at( mappedCoarsePreviousTimeStepsList.at( i ).size() - 1 ); - colIndex++; - } - } - - assert( colIndex == nbCols ); - - // Remove dependent columns of DeltaC and DeltaF - - nbRemoveCols = 1; - - while ( nbRemoveCols > 0 ) - { - nbRemoveCols = 0; - - if ( DeltaF.cols() == 1 ) - break; - - // Calculate singular value decomposition with Eigen - - Eigen::JacobiSVD svd_F( DeltaF, Eigen::ComputeThinU | Eigen::ComputeThinV ); - - vector singularValues = svd_F.singularValues(); - - for ( int i = 0; i < singularValues.rows(); i++ ) - { - if ( std::abs( singularValues( i ) ) <= singularityLimit ) - { - // Remove the column from DeltaC and DeltaF - removeColumnFromMatrix( DeltaC, i - nbRemoveCols ); - removeColumnFromMatrix( DeltaF, i - nbRemoveCols ); - - nbRemoveCols++; - } - } - - if ( nbRemoveCols ) - Info << "Response and parameter mapping: remove " << nbRemoveCols << " columns from the DeltaC and DeltaF matrices" << endl; - } - - assert( DeltaC.cols() == DeltaF.cols() ); - assert( DeltaC.cols() >= 1 ); - - // Calculate singular value decomposition with Eigen - - Eigen::JacobiSVD svd_F( DeltaF, Eigen::ComputeThinU | Eigen::ComputeThinV ); - Eigen::JacobiSVD svd_C( DeltaC, Eigen::ComputeThinU | Eigen::ComputeThinV ); - - matrix pseudoSigma_F = svd_F.singularValues().asDiagonal(); - - for ( int i = 0; i < pseudoSigma_F.cols(); i++ ) - pseudoSigma_F( i, i ) = 1.0 / pseudoSigma_F( i, i ); - - matrix pseudoDeltaF = svd_F.matrixV() * pseudoSigma_F * svd_F.matrixU().transpose(); - - // Update the design specification yk - - matrix Tk = DeltaC * pseudoDeltaF + ( I - svd_C.matrixU() * svd_C.matrixU().transpose() ) * ( I - svd_F.matrixU() * svd_F.matrixU().transpose() ); - yk = coarseResiduals.col( k ) - Tk * (fineResiduals.col( k ) - y); - } - - // Update the fine model optimum - - Info << nl << "Response and parameter mapping: update the fine model optimum by minimizing the coarse model" << endl; - - output.resize( n ); - - if ( nbCols == 0 ) - coarseModel->optimize( yk, xk, output ); - - if ( nbCols >= 1 ) - coarseModel->optimize( yk, xk, output, B, xktilde, xk ); - - if ( !coarseModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; - - xk = output; - - output.resize( m ); - - // Fine model evaluation - - fineModel->evaluate( xk, output, R ); - fineResiduals.conservativeResize( fineResiduals.rows(), fineResiduals.cols() + 1 ); - fineResiduals.col( fineResiduals.cols() - 1 ) = R; - - iter++; - - // Parameter extraction - Info << nl << "Response and parameter mapping: parameter extraction process" << endl; - coarseModel->optimize( R, xk, xktilde ); - xktildeList.push_back( xktilde ); - xkList.push_back( xk ); - - if ( !coarseModel->allConverged() ) - Warning << "Surrogate model optimization process is not converged." << endl; - - // Coarse model evaluation - - coarseModel->evaluate( xktilde, output, R ); - coarseResiduals.conservativeResize( coarseResiduals.rows(), coarseResiduals.cols() + 1 ); - coarseResiduals.col( coarseResiduals.cols() - 1 ) = R; - - assert( coarseResiduals.cols() == fineResiduals.cols() ); - assert( static_cast( xktildeList.size() ) == coarseResiduals.cols() ); - assert( static_cast( xkList.size() ) == coarseResiduals.cols() ); - - // Check convergence criteria - if ( isConvergence( xk, xkprev, residualCriterium ) ) - break; - } - - iterationsConverged(); -} - -void ResponseParameterMapping::removeColumnFromMatrix( - matrix & A, - int col - ) -{ - for ( int j = col; j < A.cols() - 1; j++ ) - A.col( j ) = A.col( j + 1 ); - - A.conservativeResize( A.rows(), A.cols() - 1 ); -} diff --git a/src/fsi/ResponseParameterMapping.H b/src/fsi/ResponseParameterMapping.H deleted file mode 100644 index 2047be05..00000000 --- a/src/fsi/ResponseParameterMapping.H +++ /dev/null @@ -1,56 +0,0 @@ - -/* - * Author - * David Blom, TU Delft. All rights reserved. - */ - -#ifndef ResponseParameterMapping_H -#define ResponseParameterMapping_H - -#include "SpaceMapping.H" - -namespace fsi -{ - class ResponseParameterMapping : public SpaceMapping - { -public: - - using SpaceMapping::performPostProcessing; - - ResponseParameterMapping( - shared_ptr fineModel, - shared_ptr coarseModel, - int maxIter, - int nbReuse, - int reuseInformationStartingFromTimeIndex, - double singularityLimit - ); - - virtual ~ResponseParameterMapping(); - - virtual void iterationsConverged(); - - virtual void performPostProcessing( - const vector & y, - const vector & x0, - vector & xk, - bool residualCriterium - ); - - void removeColumnFromMatrix( - matrix & A, - int col - ); - - int iter; - - std::deque xktildeList; - std::deque xkList; - std::deque mappedCoarseList; - std::deque > xktildePreviousTimeStepsList; - std::deque > xkPreviousTimeStepsList; - std::deque > mappedCoarsePreviousTimeStepsList; - }; -} - -#endif diff --git a/src/fsi/SpaceMapping.C b/src/fsi/SpaceMapping.C index 80f5400e..842a20cf 100644 --- a/src/fsi/SpaceMapping.C +++ b/src/fsi/SpaceMapping.C @@ -12,6 +12,7 @@ SpaceMapping::SpaceMapping( shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit @@ -20,6 +21,7 @@ SpaceMapping::SpaceMapping( fineModel( fineModel ), coarseModel( coarseModel ), maxIter( maxIter ), + maxUsedIterations( maxUsedIterations ), nbReuse( nbReuse ), reuseInformationStartingFromTimeIndex( reuseInformationStartingFromTimeIndex ), singularityLimit( singularityLimit ), @@ -27,7 +29,9 @@ SpaceMapping::SpaceMapping( coarseResiduals(), fineResiduals(), coarseResidualsList(), - fineResidualsList() + fineResidualsList(), + coarseResidualsTimeList(), + fineResidualsTimeList() { assert( fineModel ); assert( coarseModel ); @@ -36,10 +40,35 @@ SpaceMapping::SpaceMapping( assert( reuseInformationStartingFromTimeIndex >= 0 ); assert( singularityLimit > 0 ); assert( singularityLimit < 1 ); + assert( maxUsedIterations > 0 ); }; void SpaceMapping::finalizeTimeStep() { + // Save input/output information for next time step + if ( nbReuse > 0 && fineResidualsList.size() >= 1 && timeIndex >= reuseInformationStartingFromTimeIndex ) + { + assert( fineResidualsTimeList.size() == coarseResidualsTimeList.size() ); + fineResidualsTimeList.push_front( fineResidualsList ); + coarseResidualsTimeList.push_front( coarseResidualsList ); + } + + // Remove the last items from the residual list and solutions list + // in order to ensure that at maximum nbReuse time steps + // are included. + while ( static_cast( fineResidualsTimeList.size() ) > nbReuse ) + { + fineResidualsTimeList.pop_back(); + coarseResidualsTimeList.pop_back(); + } + + fineResidualsList.clear(); + coarseResidualsList.clear(); + + assert( fineResidualsTimeList.size() == coarseResidualsTimeList.size() ); + assert( coarseResiduals.size() == fineResiduals.size() ); + assert( static_cast( fineResidualsTimeList.size() ) <= nbReuse ); + timeIndex++; } @@ -57,29 +86,19 @@ bool SpaceMapping::isConvergence( void SpaceMapping::iterationsConverged() { - // Save input/output information for next time step - if ( nbReuse > 0 && fineResiduals.cols() >= 2 && timeIndex >= reuseInformationStartingFromTimeIndex ) + // Save input/output information for next solve + if ( fineResiduals.size() >= 2 ) { - assert( coarseResiduals.cols() == fineResiduals.cols() ); + assert( fineResiduals.size() == coarseResiduals.size() ); - coarseResidualsList.push_front( coarseResiduals ); fineResidualsList.push_front( fineResiduals ); + coarseResidualsList.push_front( coarseResiduals ); } - // Remove the last item from the residual list and solutions - // list in order to ensure that at maximum nbReuse time steps - // are included. - while ( static_cast( fineResidualsList.size() ) > nbReuse ) - { - coarseResidualsList.pop_back(); - fineResidualsList.pop_back(); - } - - coarseResiduals.resize( coarseResiduals.rows(), 0 ); - fineResiduals.resize( fineResiduals.rows(), 0 ); + fineResiduals.clear(); + coarseResiduals.clear(); - assert( coarseResidualsList.size() == fineResidualsList.size() ); - assert( static_cast( coarseResidualsList.size() ) <= nbReuse ); + assert( fineResidualsList.size() == coarseResidualsList.size() ); } void SpaceMapping::performPostProcessing( diff --git a/src/fsi/SpaceMapping.H b/src/fsi/SpaceMapping.H index e9ec1449..99cfa939 100644 --- a/src/fsi/SpaceMapping.H +++ b/src/fsi/SpaceMapping.H @@ -13,6 +13,7 @@ #include "fvCFD.H" using namespace fsi; +using std::deque; class SpaceMapping { @@ -22,6 +23,7 @@ public: shared_ptr fineModel, shared_ptr coarseModel, int maxIter, + int maxUsedIterations, int nbReuse, int reuseInformationStartingFromTimeIndex, double singularityLimit @@ -61,15 +63,18 @@ public: shared_ptr coarseModel; int maxIter; + int maxUsedIterations; int nbReuse; int reuseInformationStartingFromTimeIndex; double singularityLimit; int timeIndex; - matrix coarseResiduals; - matrix fineResiduals; - std::deque coarseResidualsList; - std::deque fineResidualsList; + deque coarseResiduals; + deque fineResiduals; + deque > coarseResidualsList; + deque > fineResidualsList; + deque > > coarseResidualsTimeList; + deque > > fineResidualsTimeList; }; #endif diff --git a/src/tests/app/Make/files b/src/tests/app/Make/files index 1a69fdcd..6de6be57 100644 --- a/src/tests/app/Make/files +++ b/src/tests/app/Make/files @@ -59,5 +59,6 @@ test_linearizedfluidsolver.C test_manifoldmappinglinearizedfsisolver.C test_rbfcoarsening.C test_multilevelspacemapping.C +test_multilevelmliqnilssolver.C EXE = $(FOAM_USER_APPBIN)/tests diff --git a/src/tests/app/ResponseParameterMapping.C b/src/tests/app/ResponseParameterMapping.C deleted file mode 120000 index 9aada83e..00000000 --- a/src/tests/app/ResponseParameterMapping.C +++ /dev/null @@ -1 +0,0 @@ -../../fsi/ResponseParameterMapping.C \ No newline at end of file diff --git a/src/tests/app/ResponseParameterMapping.H b/src/tests/app/ResponseParameterMapping.H deleted file mode 120000 index ecf65d90..00000000 --- a/src/tests/app/ResponseParameterMapping.H +++ /dev/null @@ -1 +0,0 @@ -../../fsi/ResponseParameterMapping.H \ No newline at end of file diff --git a/src/tests/app/test_aggressivespacemapping.C b/src/tests/app/test_aggressivespacemapping.C index 1b2bf61f..d66e8dcb 100644 --- a/src/tests/app/test_aggressivespacemapping.C +++ b/src/tests/app/test_aggressivespacemapping.C @@ -91,9 +91,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); @@ -173,9 +170,14 @@ protected: shared_ptr coarseModel( new ImplicitMultiLevelFsiSolver( multiLevelFsiSolver, postProcessing ) ); - // Create manifold mapping object + // Create aggressive space mapping object + + maxUsedIterations = couplingGridSize; + + if ( parallel ) + maxUsedIterations *= 2; - shared_ptr aggressiveSpaceMapping( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + shared_ptr aggressiveSpaceMapping( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, aggressiveSpaceMapping ); @@ -247,8 +249,8 @@ TEST_P( AggressiveSpaceMappingSolverParametrizedTest, reuse ) if ( !parallel && couplingGridSize == 20 && nbReuse == 4 && extrapolation == 0 && minIter == 3 ) { - ASSERT_EQ( solver->fineModel->fsi->nbIter, 567 ); - ASSERT_EQ( solver->coarseModel->fsi->nbIter, 3298 ); + ASSERT_EQ( solver->fineModel->fsi->nbIter, 543 ); + ASSERT_EQ( solver->coarseModel->fsi->nbIter, 3101 ); } } @@ -400,7 +402,7 @@ protected: // Create manifold mapping object - shared_ptr aggressiveSpaceMapping( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + shared_ptr aggressiveSpaceMapping( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, aggressiveSpaceMapping ); diff --git a/src/tests/app/test_andersonpostprocessing.C b/src/tests/app/test_andersonpostprocessing.C index e7e514c3..5a8fd57b 100644 --- a/src/tests/app/test_andersonpostprocessing.C +++ b/src/tests/app/test_andersonpostprocessing.C @@ -79,9 +79,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); diff --git a/src/tests/app/test_asmils.C b/src/tests/app/test_asmils.C index e4cdcee1..5d1669c9 100644 --- a/src/tests/app/test_asmils.C +++ b/src/tests/app/test_asmils.C @@ -49,7 +49,7 @@ protected: double tol = 1.0e-5; int maxIter = 500; double initialRelaxation = 1.0e-3; - double singularityLimit = 1.0e-11; + double singularityLimit = 1.0e-12; int reuseInformationStartingFromTimeIndex = 0; bool scaling = false; double beta = 1; @@ -91,9 +91,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); @@ -173,9 +170,14 @@ protected: shared_ptr coarseModel( new ImplicitMultiLevelFsiSolver( multiLevelFsiSolver, postProcessing ) ); - // Create manifold mapping object + // Create ASMILS object + + maxUsedIterations = couplingGridSize; + + if ( parallel ) + maxUsedIterations *= 2; - shared_ptr aggressiveSpaceMapping( new ASMILS( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + shared_ptr aggressiveSpaceMapping( new ASMILS( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, aggressiveSpaceMapping ); @@ -379,7 +381,7 @@ protected: // Create manifold mapping object - shared_ptr aggressiveSpaceMapping( new ASMILS( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + shared_ptr aggressiveSpaceMapping( new ASMILS( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, aggressiveSpaceMapping ); diff --git a/src/tests/app/test_manifoldmapping.C b/src/tests/app/test_manifoldmapping.C index 61e6bb6e..8b226232 100644 --- a/src/tests/app/test_manifoldmapping.C +++ b/src/tests/app/test_manifoldmapping.C @@ -426,13 +426,9 @@ TEST( ManifoldMapping, specification1 ) // Settings int maxIter = 200; - double singularityLimit = 1.0e-15; - int nbReuse = 0; - int reuseInformationStartingFromTimeIndex = 0; - bool updateJacobian = false; // Create manifold mapping object - ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ); + ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter ); mm.performPostProcessing( y, x0, xk ); @@ -458,13 +454,9 @@ TEST( ManifoldMapping, specification2 ) // Settings int maxIter = 500; - double singularityLimit = 1.0e-15; - int nbReuse = 0; - int reuseInformationStartingFromTimeIndex = 0; - bool updateJacobian = false; // Create manifold mapping object - ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ); + ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter ); mm.performPostProcessing( y, x0, xk ); @@ -490,13 +482,9 @@ TEST( ManifoldMapping, specification3 ) // Settings int maxIter = 200; - double singularityLimit = 1.0e-15; - int nbReuse = 0; - int reuseInformationStartingFromTimeIndex = 0; - bool updateJacobian = false; // Create manifold mapping object - ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ); + ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter ); mm.performPostProcessing( y, x0, xk ); @@ -522,13 +510,9 @@ TEST( ManifoldMapping, specification4 ) // Settings int maxIter = 2000; - double singularityLimit = 1.0e-15; - int nbReuse = 0; - int reuseInformationStartingFromTimeIndex = 0; - bool updateJacobian = false; // Create manifold mapping object - ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ); + ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter ); mm.performPostProcessing( y, x0, xk ); @@ -554,13 +538,9 @@ TEST( ManifoldMapping, specification5 ) // Settings int maxIter = 200; - double singularityLimit = 1.0e-15; - int nbReuse = 0; - int reuseInformationStartingFromTimeIndex = 0; - bool updateJacobian = false; // Create manifold mapping object - ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ); + ManifoldMapping mm = ManifoldMapping( fineModel, coarseModel, maxIter ); mm.performPostProcessing( y, x0, xk ); diff --git a/src/tests/app/test_manifoldmappinglinearizedfsisolver.C b/src/tests/app/test_manifoldmappinglinearizedfsisolver.C index a6d0f353..36f27d38 100644 --- a/src/tests/app/test_manifoldmappinglinearizedfsisolver.C +++ b/src/tests/app/test_manifoldmappinglinearizedfsisolver.C @@ -55,7 +55,7 @@ protected: double tol = 1.0e-5; int maxIter = 50; double initialRelaxation = 1.0e-3; - double singularityLimit = 1.0e-11; + double singularityLimit = 1.0e-12; int reuseInformationStartingFromTimeIndex = 0; bool scaling = false; double beta = 1; @@ -186,7 +186,7 @@ protected: // Create manifold mapping object - shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, manifoldMapping ); diff --git a/src/tests/app/test_manifoldmappingsolver.C b/src/tests/app/test_manifoldmappingsolver.C index 64cc4758..41d8d6cb 100644 --- a/src/tests/app/test_manifoldmappingsolver.C +++ b/src/tests/app/test_manifoldmappingsolver.C @@ -96,9 +96,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); @@ -189,7 +186,12 @@ protected: // Create manifold mapping object - shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + maxUsedIterations = couplingGridSize; + + if ( parallel ) + maxUsedIterations *= 2; + + shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, manifoldMapping ); @@ -327,9 +329,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); @@ -408,7 +407,7 @@ protected: // Create manifold mapping object - shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, manifoldMapping ); diff --git a/src/tests/app/test_mliqnilssolver.C b/src/tests/app/test_mliqnilssolver.C index 29c09d0f..b97627a1 100644 --- a/src/tests/app/test_mliqnilssolver.C +++ b/src/tests/app/test_mliqnilssolver.C @@ -97,9 +97,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); @@ -225,6 +222,14 @@ TEST_P( MLIQNILSSolverParametrizedTest, object ) ASSERT_TRUE( true ); } +TEST_P( MLIQNILSSolverParametrizedTest, solveTimeStep ) +{ + solver->solveTimeStep(); + + ASSERT_TRUE( solver->fineModel->fsi->allConverged ); + ASSERT_TRUE( solver->models->at( 0 )->fsi->allConverged ); +} + TEST_P( MLIQNILSSolverParametrizedTest, run ) { solver->run(); @@ -344,9 +349,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); diff --git a/src/tests/app/test_multilevelaggressivespacemappingsolver.C b/src/tests/app/test_multilevelaggressivespacemappingsolver.C index be9c32a7..6b2aa266 100644 --- a/src/tests/app/test_multilevelaggressivespacemappingsolver.C +++ b/src/tests/app/test_multilevelaggressivespacemappingsolver.C @@ -193,7 +193,7 @@ protected: // Create manifold mapping object - shared_ptr aggressiveSpaceMapping( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + shared_ptr aggressiveSpaceMapping( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver shared_ptr aggressiveSpaceMappingSolver( new SpaceMappingSolver( fineModel, coarseModel, aggressiveSpaceMapping ) ); @@ -255,7 +255,7 @@ protected: // Create manifold mapping object - aggressiveSpaceMapping = shared_ptr ( new AggressiveSpaceMapping( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + aggressiveSpaceMapping = shared_ptr ( new AggressiveSpaceMapping( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver aggressiveSpaceMappingSolver = shared_ptr ( new SpaceMappingSolver( fineModel, coarseModelSolver, aggressiveSpaceMapping ) ); @@ -314,7 +314,6 @@ TEST_P( MultiLevelAggressiveSpaceMappingSolverParametrizedTest, monolithic ) else ASSERT_FALSE( solver->models->at( solver->models->size() - 1 )->fsi->fluid->isRunning() ); - ASSERT_TRUE( solver->models->at( 0 )->fsi->allConverged ); ASSERT_TRUE( solver->models->at( 1 )->fsi->allConverged ); ASSERT_TRUE( solver->models->at( 2 )->fsi->allConverged ); ASSERT_NEAR( solver->models->at( solver->models->size() - 1 )->fsi->fluid->data.norm(), monolithicSolver->pn.norm(), tol ); @@ -354,15 +353,15 @@ TEST_P( MultiLevelAggressiveSpaceMappingSolverParametrizedTest, timeStep ) if ( couplingGridSize == 50 && parallel ) { - ASSERT_EQ( solver->models->at( 0 )->fsi->nbIter, 506 ); - ASSERT_EQ( solver->models->at( 1 )->fsi->nbIter, 59 ); + ASSERT_EQ( solver->models->at( 0 )->fsi->nbIter, 446 ); + ASSERT_EQ( solver->models->at( 1 )->fsi->nbIter, 53 ); ASSERT_EQ( solver->models->at( 2 )->fsi->nbIter, 6 ); } if ( couplingGridSize == 50 && !parallel ) { - ASSERT_EQ( solver->models->at( 0 )->fsi->nbIter, 185 ); - ASSERT_EQ( solver->models->at( 1 )->fsi->nbIter, 31 ); + ASSERT_EQ( solver->models->at( 0 )->fsi->nbIter, 142 ); + ASSERT_EQ( solver->models->at( 1 )->fsi->nbIter, 25 ); ASSERT_EQ( solver->models->at( 2 )->fsi->nbIter, 4 ); } } diff --git a/src/tests/app/test_multilevelasmils.C b/src/tests/app/test_multilevelasmils.C index b7fe100e..2578b790 100644 --- a/src/tests/app/test_multilevelasmils.C +++ b/src/tests/app/test_multilevelasmils.C @@ -193,7 +193,7 @@ protected: // Create manifold mapping object - shared_ptr aggressiveSpaceMapping( new ASMILS( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + shared_ptr aggressiveSpaceMapping( new ASMILS( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver shared_ptr aggressiveSpaceMappingSolver( new SpaceMappingSolver( fineModel, coarseModel, aggressiveSpaceMapping ) ); @@ -255,7 +255,7 @@ protected: // Create manifold mapping object - aggressiveSpaceMapping = shared_ptr ( new ASMILS( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + aggressiveSpaceMapping = shared_ptr ( new ASMILS( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver aggressiveSpaceMappingSolver = shared_ptr ( new SpaceMappingSolver( fineModel, coarseModelSolver, aggressiveSpaceMapping ) ); diff --git a/src/tests/app/test_multilevelmanifoldmappingsolver.C b/src/tests/app/test_multilevelmanifoldmappingsolver.C index b308ddd0..e70bc2c6 100644 --- a/src/tests/app/test_multilevelmanifoldmappingsolver.C +++ b/src/tests/app/test_multilevelmanifoldmappingsolver.C @@ -102,9 +102,6 @@ protected: shared_ptr rbfInterpToCouplingMesh; shared_ptr rbfInterpToMesh; - - - rbfFunction = shared_ptr( new TPSFunction() ); rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); @@ -194,7 +191,7 @@ protected: // Create manifold mapping object - shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + shared_ptr manifoldMapping( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); // Create manifold mapping solver shared_ptr manifoldMappingSolver( new SpaceMappingSolver( fineModel, coarseModel, manifoldMapping ) ); @@ -256,7 +253,7 @@ protected: // Create manifold mapping object - manifoldMapping = shared_ptr ( new ManifoldMapping( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + manifoldMapping = shared_ptr ( new ManifoldMapping( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); // Create manifold mapping solver manifoldMappingSolver = shared_ptr ( new SpaceMappingSolver( fineModel, coarseModelSolver, manifoldMapping ) ); @@ -325,7 +322,6 @@ TEST_P( MultiLevelManifoldMappingSolverParametrizedTest, monolithic ) else ASSERT_FALSE( solver->models->at( solver->models->size() - 1 )->fsi->fluid->isRunning() ); - ASSERT_TRUE( solver->models->at( 0 )->fsi->allConverged ); ASSERT_TRUE( solver->models->at( 1 )->fsi->allConverged ); ASSERT_TRUE( solver->models->at( 2 )->fsi->allConverged ); ASSERT_NEAR( solver->models->at( solver->models->size() - 1 )->fsi->fluid->data.norm(), monolithicSolver->pn.norm(), tol ); diff --git a/src/tests/app/test_multilevelmliqnilssolver.C b/src/tests/app/test_multilevelmliqnilssolver.C new file mode 100644 index 00000000..44974e38 --- /dev/null +++ b/src/tests/app/test_multilevelmliqnilssolver.C @@ -0,0 +1,272 @@ + +/* + * Author + * David Blom, TU Delft. All rights reserved. + */ + +#include "AndersonPostProcessing.H" +#include "MLIQNILSSolver.H" +#include "MinIterationConvergenceMeasure.H" +#include "MonolithicFsiSolver.H" +#include "MultiLevelFsiSolver.H" +#include "MultiLevelSpaceMappingSolver.H" +#include "RelativeConvergenceMeasure.H" +#include "SpaceMappingSolver.H" +#include "TubeFlowFluidSolver.H" +#include "TubeFlowSolidSolver.H" +#include "gtest/gtest.h" + +using namespace tubeflow; +using::testing::TestWithParam; +using::testing::Bool; +using::testing::Values; +using::testing::Combine; + +class MultiLevelMLIQNILSSolverParametrizedTest : public TestWithParam< std::tr1::tuple > +{ +protected: + + virtual void SetUp() + { + // Physical settings + double r0 = 0.2; + double a0 = M_PI * r0 * r0; + double u0 = 0.1; + double p0 = 0; + double dt = 0.1; + int N = 10; + double L = 1; + double T = 10; + double dx = L / N; + double rho = 1.225; + double E = 490; + double h = 1.0e-3; + double cmk = std::sqrt( E * h / (2 * rho * r0) ); + double c0 = std::sqrt( cmk * cmk - p0 / (2 * rho) ); + double kappa = c0 / u0; + double tau = u0 * dt / L; + + // Computational settings + double tol = 1.0e-5; + int maxIter = 500; + double initialRelaxation = 1.0e-3; + double singularityLimit = 1.0e-12; + bool scaling = false; + bool updateJacobian = false; + double beta = 1; + + // Parametrized settings + bool parallel = std::tr1::get<0>( GetParam() ); + int nbReuse = std::tr1::get<1>( GetParam() ); + int extrapolation = std::tr1::get<2>( GetParam() ); + int minIter = std::tr1::get<3>( GetParam() ); + int couplingGridSize = std::tr1::get<4>( GetParam() ); + int reuseInformationStartingFromTimeIndex = 0; + bool convergenceMeasureTraction = std::tr1::get<5>( GetParam() ); + + ASSERT_NEAR( tau, 0.01, 1.0e-13 ); + ASSERT_NEAR( kappa, 10, 1.0e-13 ); + ASSERT_TRUE( dx > 0 ); + + // Initialize variables + shared_ptr fluid; + shared_ptr solid; + shared_ptr fineModelFluid; + shared_ptr fineModelSolid; + shared_ptr multiLevelFluidSolver; + shared_ptr multiLevelSolidSolver; + shared_ptr< std::list > > convergenceMeasures; + shared_ptr multiLevelFsiSolver; + shared_ptr postProcessing; + shared_ptr< std::deque > > models; + + models = shared_ptr< std::deque > > ( new std::deque > () ); + + fineModelFluid = shared_ptr ( new TubeFlowFluidSolver( a0, u0, p0, dt, cmk, couplingGridSize, L, T, rho ) ); + fineModelSolid = shared_ptr ( new TubeFlowSolidSolver( a0, cmk, p0, rho, L, couplingGridSize ) ); + + // Coarse model + + int maxUsedIterations = N; + + if ( parallel ) + maxUsedIterations *= 2; + + fluid = shared_ptr ( new TubeFlowFluidSolver( a0, u0, p0, dt, cmk, N, L, T, rho ) ); + solid = shared_ptr ( new TubeFlowSolidSolver( a0, cmk, p0, rho, L, N ) ); + + shared_ptr rbfFunction; + shared_ptr rbfInterpolator; + shared_ptr rbfInterpToCouplingMesh; + shared_ptr rbfInterpToMesh; + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + multiLevelFluidSolver = shared_ptr ( new MultiLevelSolver( fluid, fineModelFluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 0, 0 ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + multiLevelSolidSolver = shared_ptr ( new MultiLevelSolver( solid, fineModelSolid, rbfInterpToCouplingMesh, rbfInterpToMesh, 1, 0 ) ); + + // Convergence measures + convergenceMeasures = shared_ptr > >( new std::list > ); + + convergenceMeasures->push_back( shared_ptr ( new MinIterationConvergenceMeasure( 0, 5 ) ) ); + convergenceMeasures->push_back( shared_ptr ( new RelativeConvergenceMeasure( 0, 1.0e-3 * tol ) ) ); + + if ( parallel || convergenceMeasureTraction ) + convergenceMeasures->push_back( std::shared_ptr ( new RelativeConvergenceMeasure( 1, 1.0e-3 * tol ) ) ); + + multiLevelFsiSolver = shared_ptr ( new MultiLevelFsiSolver( multiLevelFluidSolver, multiLevelSolidSolver, convergenceMeasures, parallel, extrapolation ) ); + + postProcessing = shared_ptr ( new AndersonPostProcessing( multiLevelFsiSolver, maxIter, initialRelaxation, maxUsedIterations, nbReuse, singularityLimit, reuseInformationStartingFromTimeIndex, scaling, beta, updateJacobian ) ); + + shared_ptr coarseModel( new ImplicitMultiLevelFsiSolver( multiLevelFsiSolver, postProcessing ) ); + + fluid.reset(); + solid.reset(); + multiLevelFluidSolver.reset(); + multiLevelSolidSolver.reset(); + convergenceMeasures.reset(); + multiLevelFsiSolver.reset(); + postProcessing.reset(); + + // Fine model + + maxUsedIterations = 2 * N; + + if ( parallel ) + maxUsedIterations *= 2; + + fluid = shared_ptr ( new TubeFlowFluidSolver( a0, u0, p0, dt, cmk, 2 * N, L, T, rho ) ); + solid = shared_ptr ( new TubeFlowSolidSolver( a0, cmk, p0, rho, L, 2 * N ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + multiLevelFluidSolver = shared_ptr ( new MultiLevelSolver( fluid, fineModelFluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 0, 1 ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + multiLevelSolidSolver = shared_ptr ( new MultiLevelSolver( solid, fineModelSolid, rbfInterpToCouplingMesh, rbfInterpToMesh, 1, 1 ) ); + + // Convergence measures + convergenceMeasures = shared_ptr > >( new std::list > ); + + convergenceMeasures->push_back( shared_ptr ( new MinIterationConvergenceMeasure( 0, minIter ) ) ); + convergenceMeasures->push_back( shared_ptr ( new RelativeConvergenceMeasure( 0, 1.0e-2 * tol ) ) ); + + if ( parallel || convergenceMeasureTraction ) + convergenceMeasures->push_back( std::shared_ptr ( new RelativeConvergenceMeasure( 1, 1.0e-2 * tol ) ) ); + + multiLevelFsiSolver = shared_ptr ( new MultiLevelFsiSolver( multiLevelFluidSolver, multiLevelSolidSolver, convergenceMeasures, parallel, extrapolation ) ); + postProcessing = shared_ptr ( new AndersonPostProcessing( multiLevelFsiSolver, maxIter, initialRelaxation, maxUsedIterations, nbReuse, singularityLimit, reuseInformationStartingFromTimeIndex, scaling, beta, updateJacobian ) ); + + shared_ptr fineModel( new ImplicitMultiLevelFsiSolver( multiLevelFsiSolver, postProcessing ) ); + + models->push_back( coarseModel ); + models->push_back( fineModel ); + + fluid.reset(); + solid.reset(); + multiLevelFluidSolver.reset(); + multiLevelSolidSolver.reset(); + convergenceMeasures.reset(); + multiLevelFsiSolver.reset(); + postProcessing.reset(); + + // Fine model 2 + + maxUsedIterations = couplingGridSize; + + if ( parallel ) + maxUsedIterations *= 2; + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + multiLevelFluidSolver = shared_ptr ( new MultiLevelSolver( fineModelFluid, fineModelFluid, rbfInterpToCouplingMesh, rbfInterpToMesh, 0, 2 ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToCouplingMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + rbfFunction = shared_ptr( new TPSFunction() ); + rbfInterpolator = shared_ptr( new RBFInterpolation( rbfFunction ) ); + rbfInterpToMesh = shared_ptr ( new RBFCoarsening( rbfInterpolator ) ); + + multiLevelSolidSolver = shared_ptr ( new MultiLevelSolver( fineModelSolid, fineModelSolid, rbfInterpToCouplingMesh, rbfInterpToMesh, 1, 2 ) ); + + // Convergence measures + convergenceMeasures = shared_ptr > >( new std::list > ); + + convergenceMeasures->push_back( shared_ptr ( new MinIterationConvergenceMeasure( 0, 1 ) ) ); + convergenceMeasures->push_back( shared_ptr ( new RelativeConvergenceMeasure( 0, tol ) ) ); + + if ( parallel || convergenceMeasureTraction ) + convergenceMeasures->push_back( std::shared_ptr ( new RelativeConvergenceMeasure( 1, tol ) ) ); + + multiLevelFsiSolver = shared_ptr ( new MultiLevelFsiSolver( multiLevelFluidSolver, multiLevelSolidSolver, convergenceMeasures, parallel, extrapolation ) ); + postProcessing = shared_ptr ( new AndersonPostProcessing( multiLevelFsiSolver, maxIter, initialRelaxation, maxUsedIterations, nbReuse, singularityLimit, reuseInformationStartingFromTimeIndex, scaling, beta, updateJacobian ) ); + + fineModel = shared_ptr ( new ImplicitMultiLevelFsiSolver( multiLevelFsiSolver, postProcessing ) ); + + models->push_back( fineModel ); + + // Create multi level manifold mapping solver + + solver = new MLIQNILSSolver( models, true ); + + // Monolithic solver + monolithicSolver = new MonolithicFsiSolver( a0, u0, p0, dt, cmk, couplingGridSize, L, T, rho ); + } + + virtual void TearDown() + { + delete solver; + delete monolithicSolver; + } + + MonolithicFsiSolver * monolithicSolver; + MLIQNILSSolver * solver; +}; + +INSTANTIATE_TEST_CASE_P( testParameters, MultiLevelMLIQNILSSolverParametrizedTest, ::testing::Combine( Bool(), Values( 0 ), Values( 0 ), Values( 3 ), Values( 20, 40 ), Bool() ) ); + +TEST_P( MultiLevelMLIQNILSSolverParametrizedTest, object ) +{ + ASSERT_TRUE( true ); +} + +TEST_P( MultiLevelMLIQNILSSolverParametrizedTest, timeStep ) +{ + solver->solveTimeStep(); +} diff --git a/src/tests/app/test_multileveloutputspacemappingsolver.C b/src/tests/app/test_multileveloutputspacemappingsolver.C index d88a0fd3..cbef0058 100644 --- a/src/tests/app/test_multileveloutputspacemappingsolver.C +++ b/src/tests/app/test_multileveloutputspacemappingsolver.C @@ -195,7 +195,7 @@ protected: // Create manifold mapping object - shared_ptr outputSpaceMapping( new OutputSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); + shared_ptr outputSpaceMapping( new OutputSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); // Create manifold mapping solver shared_ptr outputSpaceMappingSolver( new SpaceMappingSolver( fineModel, coarseModel, outputSpaceMapping ) ); @@ -257,7 +257,7 @@ protected: // Create manifold mapping object - outputSpaceMapping = shared_ptr ( new OutputSpaceMapping( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); + outputSpaceMapping = shared_ptr ( new OutputSpaceMapping( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); // Create manifold mapping solver outputSpaceMappingSolver = shared_ptr ( new SpaceMappingSolver( fineModel, coarseModelSolver, outputSpaceMapping ) ); diff --git a/src/tests/app/test_multilevelspacemapping.C b/src/tests/app/test_multilevelspacemapping.C index 64c943f1..79bb0d27 100644 --- a/src/tests/app/test_multilevelspacemapping.C +++ b/src/tests/app/test_multilevelspacemapping.C @@ -23,7 +23,7 @@ using::testing::Bool; using::testing::Values; using::testing::Combine; -class MultiLevelSpaceMappingSolverParametrizedTest : public TestWithParam< std::tr1::tuple > +class MultiLevelSpaceMappingSolverParametrizedTest : public TestWithParam< std::tr1::tuple > { protected: @@ -68,8 +68,8 @@ protected: int reuseInformationStartingFromTimeIndex = 0; bool convergenceMeasureTraction = std::tr1::get<5>( GetParam() ); int spaceMappingAlgorithm = std::tr1::get<6>( GetParam() ); - bool coarsening = std::tr1::get<7>( GetParam() ); - bool liveSelection = std::tr1::get<8>( GetParam() ); + bool coarsening = false; + bool liveSelection = false; assert( spaceMappingAlgorithm > -1 && spaceMappingAlgorithm < 3 ); @@ -203,13 +203,13 @@ protected: shared_ptr spaceMapping; if ( spaceMappingAlgorithm == 0 ) - spaceMapping = shared_ptr( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + spaceMapping = shared_ptr( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); if ( spaceMappingAlgorithm == 1 ) - spaceMapping = shared_ptr( new OutputSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + spaceMapping = shared_ptr( new OutputSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); if ( spaceMappingAlgorithm == 2 ) - spaceMapping = shared_ptr( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + spaceMapping = shared_ptr( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); assert( spaceMapping ); @@ -274,13 +274,13 @@ protected: // Create space mapping object if ( spaceMappingAlgorithm == 0 ) - spaceMapping = shared_ptr( new ManifoldMapping( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + spaceMapping = shared_ptr( new ManifoldMapping( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); if ( spaceMappingAlgorithm == 1 ) - spaceMapping = shared_ptr( new OutputSpaceMapping( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + spaceMapping = shared_ptr( new OutputSpaceMapping( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); if ( spaceMappingAlgorithm == 2 ) - spaceMapping = shared_ptr( new AggressiveSpaceMapping( fineModel, coarseModelSolver, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + spaceMapping = shared_ptr( new AggressiveSpaceMapping( fineModel, coarseModelSolver, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); // Create manifold mapping solver spaceMappingSolver = shared_ptr ( new SpaceMappingSolver( fineModel, coarseModelSolver, spaceMapping ) ); @@ -301,7 +301,7 @@ protected: shared_ptr solver; }; -INSTANTIATE_TEST_CASE_P( testParameters, MultiLevelSpaceMappingSolverParametrizedTest, ::testing::Combine( Bool(), Values( 0, 4 ), Values( 0, 2 ), Values( 3 ), Values( 20, 40 ), Bool(), Values( 0, 1, 2 ), Bool(), Bool() ) ); +INSTANTIATE_TEST_CASE_P( testParameters, MultiLevelSpaceMappingSolverParametrizedTest, ::testing::Combine( Bool(), Values( 0, 4 ), Values( 0, 2 ), Values( 3 ), Values( 20, 40 ), Bool(), Values( 0, 1, 2 ) ) ); TEST_P( MultiLevelSpaceMappingSolverParametrizedTest, run ) { diff --git a/src/tests/app/test_outputspacemapping.C b/src/tests/app/test_outputspacemapping.C index 5f262f08..95f0f27f 100644 --- a/src/tests/app/test_outputspacemapping.C +++ b/src/tests/app/test_outputspacemapping.C @@ -52,7 +52,7 @@ protected: int maxIter = 500; double initialRelaxation = 1.0e-3; int reuseInformationStartingFromTimeIndex = 0; - double singularityLimit = 1.0e-10; + double singularityLimit = 1.0e-13; bool scaling = false; double beta = 1; bool updateJacobian = false; @@ -175,7 +175,7 @@ protected: // Create manifold mapping object - shared_ptr outputSpaceMapping( new OutputSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); + shared_ptr outputSpaceMapping( new OutputSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, outputSpaceMapping ); @@ -417,7 +417,7 @@ protected: // Create manifold mapping object - shared_ptr outputSpaceMapping( new OutputSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); + shared_ptr outputSpaceMapping( new OutputSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, order ) ); // Create manifold mapping solver solver = new SpaceMappingSolver( fineModel, coarseModel, outputSpaceMapping ); diff --git a/src/tests/tubeflow/tubeflow.C b/src/tests/tubeflow/tubeflow.C index 2e8c723f..5e75bd4f 100644 --- a/src/tests/tubeflow/tubeflow.C +++ b/src/tests/tubeflow/tubeflow.C @@ -333,16 +333,16 @@ int main( shared_ptr spaceMapping; if ( fsiSolver == "MM" ) - spaceMapping = shared_ptr ( new ManifoldMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); + spaceMapping = shared_ptr ( new ManifoldMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, updateJacobian ) ); if ( fsiSolver == "OSM" ) - spaceMapping = shared_ptr ( new OutputSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, iOrder ) ); + spaceMapping = shared_ptr ( new OutputSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit, iOrder ) ); if ( fsiSolver == "ASM" ) - spaceMapping = shared_ptr ( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + spaceMapping = shared_ptr ( new AggressiveSpaceMapping( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); if ( fsiSolver == "ASM-ILS" ) - spaceMapping = shared_ptr ( new ASMILS( fineModel, coarseModel, maxIter, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); + spaceMapping = shared_ptr ( new ASMILS( fineModel, coarseModel, maxIter, maxUsedIterations, nbReuse, reuseInformationStartingFromTimeIndex, singularityLimit ) ); shared_ptr spaceMappingSolver( new SpaceMappingSolver( fineModel, coarseModel, spaceMapping ) );