Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
NavierStokesAssFE_def.hpp
1#ifndef NAVIERSTOKESASSFE_def_hpp
2#define NAVIERSTOKESASSFE_def_hpp
3#include "NavierStokesAssFE_decl.hpp"
4
13
14/*void sxOne2D(double* x, double* res, double t, double* parameter){
15
16 res[0] = 1.;
17 res[1] = 0.;
18 return;
19}
20
21void syOne2D(double* x, double* res, double t, double* parameter){
22
23 res[0] = 0.;
24 res[1] = 1.;
25
26 return;
27}
28void sDummyFunc(double* x, double* res, double t, double* parameter){
29
30 return;
31}
32
33double OneFunction(double* x, int* parameter)
34{
35 return 1.0;
36}*/
37
38namespace FEDD {
39
40
41
42template<class SC,class LO,class GO,class NO>
43NavierStokesAssFE<SC,LO,GO,NO>::NavierStokesAssFE( const DomainConstPtr_Type &domainVelocity, std::string FETypeVelocity, const DomainConstPtr_Type &domainPressure, std::string FETypePressure, ParameterListPtr_Type parameterList ):
44NonLinearProblem<SC,LO,GO,NO>( parameterList, domainVelocity->getComm() ),
45A_(),
46pressureIDsLoc(new vec_int_Type(2)),
47u_rep_(),
48p_rep_(),
49viscosity_element_() // Added here also viscosity field
50{
51
52 this->nonLinearTolerance_ = this->parameterList_->sublist("Parameter").get("relNonLinTol",1.0e-6);
53 this->initNOXParameters();
54
55 this->addVariable( domainVelocity , FETypeVelocity , "u" , domainVelocity->getDimension());
56 this->addVariable( domainPressure , FETypePressure , "p" , 1);
57 this->dim_ = this->getDomain(0)->getDimension();
58
59 u_rep_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated() ) );
60
61 p_rep_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(1)->getMapRepeated() ) );
62
63 if (parameterList->sublist("Parameter").get("Calculate Coefficients",false)) {
64 vec2D_dbl_ptr_Type vectmpPointsPressure = domainPressure->getPointsUnique();
65 vec2D_dbl_Type::iterator it;
66 int front = -1;
67 int back = -1;
68 if (domainPressure->getDimension() == 2) {
69 it = find_if(vectmpPointsPressure->begin(), vectmpPointsPressure->end(),
70 [&] (const std::vector<double>& a){
71 if (a.at(0) >= 0.15-1.e-12 && a.at(0) <= 0.15+1.e-12
72 && a.at(1) >= 0.2-1.e-12 && a.at(1) <= 0.2+1.e-12) {
73 return true;
74 }
75 else {
76 return false;
77 }
78 });
79
80 if (it != vectmpPointsPressure->end()) {
81 front = distance(vectmpPointsPressure->begin(),it);
82 }
83 it = find_if(vectmpPointsPressure->begin(), vectmpPointsPressure->end(),
84 [&] (const std::vector<double>& a){
85 if (a.at(0) >= 0.25-1.e-12 && a.at(0) <= 0.25+1.e-12
86 && a.at(1) >= 0.2-1.e-12 && a.at(1) <= 0.2+1.e-12) {
87 return true;
88 }
89 else {
90 return false;
91 }
92 });
93
94 if (it != vectmpPointsPressure->end()) {
95 back = distance(vectmpPointsPressure->begin(),it);
96 }
97 pressureIDsLoc->at(0) = front;
98 pressureIDsLoc->at(1) = back;
99 }
100 else if(domainPressure->getDimension() == 3){
101#ifdef ASSERTS_WARNINGS
102 MYASSERT(false,"Not implemented to calc coefficients in 3D!");
103#endif
104 }
105 }
106
107}
108
109template<class SC,class LO,class GO,class NO>
110void NavierStokesAssFE<SC,LO,GO,NO>::info(){
111 this->infoProblem();
112 this->infoNonlinProblem();
113}
114
115template<class SC,class LO,class GO,class NO>
116void NavierStokesAssFE<SC,LO,GO,NO>::assemble( std::string type ) const{
117
118 if (type=="") {
119 if (this->verbose_)
120 std::cout << "-- Assembly Navier-Stokes ... " << std::endl;
121
122 assembleConstantMatrices();
123
124 if (this->verbose_)
125 std::cout << "done -- " << std::endl;
126 }
127 else
128 reAssemble( type );
129
130};
131
132template<class SC,class LO,class GO,class NO>
133void NavierStokesAssFE<SC,LO,GO,NO>::assembleConstantMatrices() const{
134
135 if (this->verbose_)
136 std::cout << "-- Assembly constant matrices Navier-Stokes ... " << std::flush;
137
138 double viscosity = this->parameterList_->sublist("Parameter").get("Viscosity",1.);
139 double density = this->parameterList_->sublist("Parameter").get("Density",1.);
140
141 // Egal welcher Wert, da OneFunction nicht von parameter abhaengt
142 int* dummy;
143
144 A_.reset(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
145
146 MapConstPtr_Type pressureMap;
147 if ( this->getDomain(1)->getFEType() == "P0" )
148 pressureMap = this->getDomain(1)->getElementMap();
149 else
150 pressureMap = this->getDomain(1)->getMapUnique();
151
152 if (this->system_.is_null())
153 this->system_.reset(new BlockMatrix_Type(2));
154
155 if (this->residualVec_.is_null())
156 this->residualVec_.reset(new BlockMultiVector_Type(2));
157
158 MatrixPtr_Type B(new Matrix_Type( pressureMap, this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
159 MatrixPtr_Type BT(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(1)->getDimension() * this->getDomain(1)->getApproxEntriesPerRow() ) );
160 MatrixPtr_Type C(new Matrix_Type( pressureMap,1));
161
162
163 this->system_->addBlock(A_,0,0);
164 this->system_->addBlock(BT,0,1);
165 this->system_->addBlock(B,1,0);
166 this->system_->addBlock(C,1,1);
167
168 this->feFactory_->assemblyNavierStokes(this->dim_, this->getDomain(0)->getFEType(), this->getDomain(1)->getFEType(), 2, this->dim_,1,u_rep_,p_rep_,this->system_,this->residualVec_,this->coeff_, this->parameterList_,false, "Jacobian", true/*call fillComplete*/);
169
170 if ( !this->getFEType(0).compare("P1") ) {
171 C.reset(new Matrix_Type( this->getDomain(1)->getMapUnique(), this->getDomain(1)->getApproxEntriesPerRow() ) );
172 this->feFactory_->assemblyBDStabilization( this->dim_, "P1", C, true);
173 C->resumeFill();
174 C->scale( -1. / ( viscosity * density ) );
175 C->fillComplete( pressureMap, pressureMap );
176
177 this->system_->addBlock( C, 1, 1 );
178 }
179
180
181
182
183#ifdef FEDD_HAVE_TEKO
184 if ( !this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic").compare("Teko") ) {
185 if (!this->parameterList_->sublist("General").get("Assemble Velocity Mass",false)) {
186 MatrixPtr_Type Mvelocity(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getApproxEntriesPerRow() ) );
187 //
188 this->feFactory_->assemblyMass( this->dim_, this->domain_FEType_vec_.at(0), "Vector", Mvelocity, true );
189 //
190 this->getPreconditionerConst()->setVelocityMassMatrix( Mvelocity );
191 if (this->verbose_)
192 std::cout << "\nVelocity mass matrix for LSC block preconditioner is assembled." << std::endl;
193 } else {
194 if (this->verbose_)
195 std::cout << "\nVelocity mass matrix for LSC block preconditioner not assembled." << std::endl;
196 }
197 }
198#endif
199 std::string precType = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
200 if ( precType == "Diagonal" || precType == "Triangular" ) {
201 MatrixPtr_Type Mpressure(new Matrix_Type( this->getDomain(1)->getMapUnique(), this->getDomain(1)->getApproxEntriesPerRow() ) );
202
203 this->feFactory_->assemblyMass( this->dim_, this->domain_FEType_vec_.at(1), "Scalar", Mpressure, true );
204 SC kinVisco = this->parameterList_->sublist("Parameter").get("Viscosity",1.);
205 Mpressure->scale(-1./kinVisco);
206 this->getPreconditionerConst()->setPressureMassMatrix( Mpressure );
207 }
208 if (this->verbose_)
209 std::cout << " Call Reassemble FixedPoint and Newton to allocate the Matrix pattern " << std::endl;
210
211 // this->reAssemble("FixedPoint");
212 this->reAssemble("Newton");
213
214 if (this->verbose_)
215 std::cout << "done -- " << std::endl;
216
217};
218
219
220template<class SC,class LO,class GO,class NO>
221void NavierStokesAssFE<SC,LO,GO,NO>::reAssemble( MatrixPtr_Type& massmatrix, std::string type ) const
222{
223
224}
225
226
227
228template<class SC,class LO,class GO,class NO>
229void NavierStokesAssFE<SC,LO,GO,NO>::reAssemble(std::string type) const {
230
231
232 if (this->verbose_)
233 std::cout << "-- Reassembly Navier-Stokes ("<< type <<") ... " << std::flush;
234
235 double density = this->parameterList_->sublist("Parameter").get("Density",1.);
236
237 MatrixPtr_Type ANW = Teuchos::rcp(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
238
239 MultiVectorConstPtr_Type u = this->solution_->getBlock(0);
240 u_rep_->importFromVector(u, true);
241
242 MultiVectorConstPtr_Type p = this->solution_->getBlock(1);
243 p_rep_->importFromVector(p, true);
244
245
246 if (type=="Rhs") {
247
248 this->system_->addBlock(ANW,0,0);
249
250 this->feFactory_->assemblyNavierStokes(this->dim_, this->getDomain(0)->getFEType(), this->getDomain(1)->getFEType(), 2, this->dim_,1,u_rep_,p_rep_,this->system_, this->residualVec_,this->coeff_,this->parameterList_, true, "FixedPoint", true); // We can also change that to "Jacobian"
251 this->feFactory_->assemblyNavierStokes(this->dim_, this->getDomain(0)->getFEType(), this->getDomain(1)->getFEType(), 2, this->dim_,1,u_rep_,p_rep_,this->system_, this->residualVec_,this->coeff_,this->parameterList_, true, "Rhs", true);
252
253 }
254 else if (type=="FixedPoint" ) {
255
256 this->system_->addBlock(ANW,0,0);
257 this->feFactory_->assemblyNavierStokes(this->dim_, this->getDomain(0)->getFEType(), this->getDomain(1)->getFEType(), 2, this->dim_,1,u_rep_,p_rep_,this->system_, this->residualVec_,this->coeff_,this->parameterList_, true, "FixedPoint", true);
258 }
259 else if(type=="Newton"){
260
261 this->system_->addBlock(ANW,0,0);
262 this->feFactory_->assemblyNavierStokes(this->dim_, this->getDomain(0)->getFEType(), this->getDomain(1)->getFEType(), 2, this->dim_,1,u_rep_,p_rep_,this->system_,this->residualVec_, this->coeff_,this->parameterList_, true,"Jacobian", true);
263
264 }
265
266 this->system_->addBlock(ANW,0,0);
267
268 if (this->verbose_)
269 std::cout << "done -- " << std::endl;
270
271
272}
273
274
275template<class SC,class LO,class GO,class NO>
276void NavierStokesAssFE<SC,LO,GO,NO>::calculateNonLinResidualVec(std::string type, double time) const{
277
278 //this->reAssemble("FixedPoint");
279 this->reAssemble("Rhs");
280
281 // We need to additionally add the residual component for the stabilization block, as it is not part of the AssembleFE routines and calculated externally here.
282 if(this->getDomain(0)->getFEType() == "P1"){
283
284 MultiVectorPtr_Type residualPressureTmp = Teuchos::rcp(new MultiVector_Type( this->getDomain(1)->getMapUnique() ));
285
286 this->system_->getBlock(1,1)->apply( *(this->solution_->getBlock(1)), *residualPressureTmp);
287
288 this->residualVec_->getBlockNonConst(1)->update(1.,*residualPressureTmp,1.);
289
290
291 }
292 // We need to account for different parameters of time discretizations here
293 // This is ok for bdf with 1.0 scaling of the system. Would be wrong for Crank-Nicolson - might be ok now for CN
294
295 /*if (this->coeff_.size() == 0)
296 this->system_->apply( *this->solution_, *this->residualVec_ );
297 else
298 this->system_->apply( *this->solution_, *this->residualVec_, this->coeff_ );*/
299
300 if (!type.compare("standard")){
301 this->residualVec_->update(-1.,*this->rhs_,1.);
302// if ( !this->sourceTerm_.is_null() )
303// this->residualVec_->update(-1.,*this->sourceTerm_,1.);
304 // this might be set again by the TimeProblem after addition of M*u
305 this->bcFactory_->setVectorMinusBC( this->residualVec_, this->solution_, time );
306
307 }
308 else if(!type.compare("reverse")){
309 this->residualVec_->update(1.,*this->rhs_,-1.); // this = -1*this + 1*rhs
310// if ( !this->sourceTerm_.is_null() )
311// this->residualVec_->update(1.,*this->sourceTerm_,1.);
312 // this might be set again by the TimeProblem after addition of M*u
313 this->bcFactory_->setBCMinusVector( this->residualVec_, this->solution_, time );
314
315 }
316
317}
318
319
320template<class SC,class LO,class GO,class NO>
321void NavierStokesAssFE<SC,LO,GO,NO>::evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
322 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs
323 ) const
324{
325 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
326 if ( !type.compare("Monolithic"))
327 evalModelImplMonolithic( inArgs, outArgs );
328 else if ( !type.compare("Teko")){
329#ifdef FEDD_HAVE_TEKO
330 evalModelImplBlock( inArgs, outArgs );
331#else
332 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
333#endif
334 }
335 else
336 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown preconditioner/solver type.");
337}
338
346template<class SC,class LO,class GO,class NO>
347void NavierStokesAssFE<SC,LO,GO,NO>::evalModelImplMonolithic(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
348 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs ) const
349{
350
351
352 using Teuchos::RCP;
353 using Teuchos::rcp;
354 using Teuchos::rcp_dynamic_cast;
355 using Teuchos::rcp_const_cast;
356 using Teuchos::ArrayView;
357 using Teuchos::Array;
358 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
359 TEUCHOS_TEST_FOR_EXCEPTION( inArgs.get_x().is_null(), std::logic_error, "inArgs.get_x() is null.");
360
361 RCP< const Thyra::VectorBase< SC > > vecThyra = inArgs.get_x();
362 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
363
364 RCP< Thyra::VectorBase< SC > > vecThyraNonConst = rcp_const_cast<Thyra::VectorBase< SC > >(vecThyra);
365
366 this->solution_->fromThyraMultiVector(vecThyraNonConst);
367
368 const RCP<Thyra::MultiVectorBase<SC> > f_out = outArgs.get_f();
369 const RCP<Thyra::LinearOpBase<SC> > W_out = outArgs.get_W_op();
370 const RCP<Thyra::PreconditionerBase<SC> > W_prec_out = outArgs.get_W_prec();
371
372 typedef Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO> tpetra_extract;
373 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraMatrix_Type;
374 typedef RCP<TpetraMatrix_Type> TpetraMatrixPtr_Type;
375 typedef RCP<const TpetraMatrix_Type> TpetraMatrixConstPtr_Type;
376
377 const bool fill_f = nonnull(f_out);
378 const bool fill_W = nonnull(W_out);
379 const bool fill_W_prec = nonnull(W_prec_out);
380
381
382 if ( fill_f || fill_W || fill_W_prec ) {
383
384 // ****************
385 // Get the underlying xpetra objects
386 // ****************
387 if (fill_f) {
388
389 this->calculateNonLinResidualVec("standard"); // Calculating residual Vector
390
391 // Changing the residualVector into a ThyraMultivector
392
393 Teuchos::RCP<Thyra::MultiVectorBase<SC> > f_thyra = this->getResidualVector()->getThyraMultiVector();
394 f_out->assign(*f_thyra);
395 }
396 TpetraMatrixPtr_Type W;
397 if (fill_W) {
398 this->reAssemble("Newton"); // ReAssembling matrices with updated u in this class
399
400 this->setBoundariesSystem(); // setting boundaries to the system
401
402 Teuchos::RCP<TpetraOp_Type> W_tpetra = tpetra_extract::getTpetraOperator(W_out);
403 Teuchos::RCP<TpetraMatrix_Type> W_tpetraMat = Teuchos::rcp_dynamic_cast<TpetraMatrix_Type>(W_tpetra);
404
405 TpetraMatrixConstPtr_Type W_systemTpetra = this->getSystem()->getMergedMatrix()->getTpetraMatrix();
406 TpetraMatrixPtr_Type W_systemTpetraNonConst = rcp_const_cast<TpetraMatrix_Type>(W_systemTpetra);
407
408 //Tpetra::CrsMatrixWrap<SC,LO,GO,NO>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<SC,LO,GO,NO>&>(*W_systemXpetraNonConst);
409 //Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>& xTpetraMat = dynamic_cast<Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>&>(*crsOp.getCrsMatrix());
410
411 Teuchos::RCP<TpetraMatrix_Type> tpetraMatTpetra = W_systemTpetraNonConst; //xTpetraMat.getTpetra_CrsMatrixNonConst();
412 W_tpetraMat->resumeFill();
413
414 for (auto i=0; i<tpetraMatTpetra->getMap()->getLocalNumElements(); i++) {
415 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type indices; //ArrayView< const LO > indices
416 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type values;
417 tpetraMatTpetra->getLocalRowView( i, indices, values);
418 W_tpetraMat->replaceLocalValues( i, indices, values);
419 }
420 W_tpetraMat->fillComplete();
421 }
422 if (fill_W_prec) {
423 this->setupPreconditioner( "Monolithic" );
424
425 // ch 26.04.19: After each setup of the preconditioner we check if we use a two-level precondtioner with multiplicative combination between the levels.
426 // If this is the case, we need to pre apply the coarse level to the residual(f_out).
427
428 std::string levelCombination = this->parameterList_->sublist("ThyraPreconditioner").sublist("Preconditioner Types").sublist("FROSch").get("Level Combination","Additive");
429 if (!levelCombination.compare("Multiplicative")) {
430 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Multiplicative Level Combination is not supported for NOX.");
431// ParameterListPtr_Type solverPList = this->getLinearSolverBuilder()->getNonconstParameterList();
432//
433// solverPList->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
434//
435// Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyra_linOp = this->getPreconditionerConst()->getThyraPrecConst()->getUnspecifiedPrecOp();
436//
437// f_out->describe(*out,Teuchos::VERB_EXTREME);
438// vecThyraNonConst->describe(*out,Teuchos::VERB_EXTREME);
439// Thyra::apply( *thyra_linOp, Thyra::NOTRANS, *f_out, vecThyraNonConst.ptr() );
440// solverPList->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",false);
441 }
442
443 }
444 }
445}
453#ifdef FEDD_HAVE_TEKO
454template<class SC,class LO,class GO,class NO>
455void NavierStokesAssFE<SC,LO,GO,NO>::evalModelImplBlock(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
456 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs ) const
457{
458
459
460 using Teuchos::RCP;
461 using Teuchos::rcp;
462 using Teuchos::rcp_dynamic_cast;
463 using Teuchos::rcp_const_cast;
464 using Teuchos::ArrayView;
465 using Teuchos::Array;
466
467 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
468 TEUCHOS_TEST_FOR_EXCEPTION( inArgs.get_x().is_null(), std::logic_error, "inArgs.get_x() is null.");
469
470 RCP< const Thyra::VectorBase< SC > > vecThyra = inArgs.get_x();
471 RCP<Teuchos::FancyOStream> out = Teuchos::VerboseObjectBase::getDefaultOStream();
472
473 RCP< Thyra::VectorBase< SC > > vecThyraNonConst = rcp_const_cast<Thyra::VectorBase< SC > >(vecThyra);
474
475 RCP< Thyra::ProductVectorBase< SC > > vecThyraBlock = rcp_dynamic_cast<Thyra::ProductVectorBase< SC > > (vecThyraNonConst);
476
477 this->solution_->getBlockNonConst(0)->fromThyraMultiVector( vecThyraBlock->getNonconstVectorBlock(0) );
478 this->solution_->getBlockNonConst(1)->fromThyraMultiVector( vecThyraBlock->getNonconstVectorBlock(1) );
479
480 const RCP<Thyra::MultiVectorBase<SC> > f_out = outArgs.get_f();
481 const RCP<Thyra::LinearOpBase<SC> > W_out = outArgs.get_W_op();
482 const RCP<Thyra::PreconditionerBase<SC> > W_prec_out = outArgs.get_W_prec();
483
484 typedef Thyra::TpetraOperatorVectorExtraction<SC,LO,GO,NO> tpetra_extract;
485 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraMatrix_Type;
486 typedef RCP<TpetraMatrix_Type> TpetraMatrixPtr_Type;
487 typedef RCP<const TpetraMatrix_Type> TpetraMatrixConstPtr_Type;
488
489 const bool fill_f = nonnull(f_out);
490 const bool fill_W = nonnull(W_out);
491 const bool fill_W_prec = nonnull(W_prec_out);
492
493 if ( fill_f || fill_W || fill_W_prec ) {
494
495 // ****************
496 // Get the underlying xpetra objects
497 // ****************
498 if (fill_f) {
499
500 this->calculateNonLinResidualVec("standard");
501
502 Teko::MultiVector f0;
503 Teko::MultiVector f1;
504 f0 = this->getResidualVector()->getBlockNonConst(0)->getThyraMultiVector();
505 f1 = this->getResidualVector()->getBlockNonConst(1)->getThyraMultiVector();
506
507 std::vector<Teko::MultiVector> f_vec; f_vec.push_back(f0); f_vec.push_back(f1);
508
509 Teko::MultiVector f = Teko::buildBlockedMultiVector(f_vec);
510
511 f_out->assign(*f);
512 }
513
514 TpetraMatrixPtr_Type W;
515 if (fill_W) {
516
517 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> TpetraCrsMatrix;
518
519 this->reAssemble("Newton");
520
521 this->setBoundariesSystem();
522
523 RCP<ThyraBlockOp_Type> W_blocks = rcp_dynamic_cast<ThyraBlockOp_Type>(W_out);
524 RCP<const ThyraOp_Type> W_block00 = W_blocks->getBlock(0,0);
525 RCP<ThyraOp_Type> W_block00NonConst = rcp_const_cast<ThyraOp_Type>( W_block00 );
526 RCP<TpetraOp_Type> W_tpetra = tpetra_extract::getTpetraOperator( W_block00NonConst );
527
528 RCP<TpetraMatrix_Type> W_tpetraMat = Teuchos::rcp_dynamic_cast<TpetraMatrix_Type>(W_tpetra);
529
530 TpetraMatrixConstPtr_Type W_matrixTpetra = this->getSystem()->getBlock(0,0)->getTpetraMatrix();
531 TpetraMatrixPtr_Type W_matrixTpetraNonConst = rcp_const_cast<TpetraMatrix_Type>(W_matrixTpetra);
532 RCP<TpetraMatrix_Type> tpetraMatTpetra = W_matrixTpetraNonConst;
533
534 W_tpetraMat->resumeFill();
535
536 for (auto i=0; i<tpetraMatTpetra->getMap()->getLocalNumElements(); i++) {
537 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::local_inds_host_view_type indices; //ArrayView< const LO > indices
538 typename Tpetra::CrsMatrix<SC,LO,GO,NO>::values_host_view_type values;
539 tpetraMatTpetra->getLocalRowView( i, indices, values);
540 W_tpetraMat->replaceLocalValues( i, indices, values);
541 }
542 W_tpetraMat->fillComplete();
543
544 }
545
546 if (fill_W_prec) {
547 if (stokesTekoPrecUsed_){
548 this->setupPreconditioner( "Teko" );
549 }
550 else
551 stokesTekoPrecUsed_ = true;
552
553 // ch 26.04.19: After each setup of the preconditioner we check if we use a two-level precondtioner with multiplicative combination between the levels.
554 // If this is the case, we need to pre apply the coarse level to the residual(f_out).
555
556 ParameterListPtr_Type tmpSubList = sublist( sublist( sublist( sublist( this->parameterList_, "Teko Parameters" ) , "Preconditioner Types" ) , "Teko" ) , "Inverse Factory Library" );
557
558 std::string levelCombination1 = tmpSubList->sublist( "FROSch-Velocity" ).get("Level Combination","Additive");
559 std::string levelCombination2 = tmpSubList->sublist( "FROSch-Pressure" ).get("Level Combination","Additive");
560
561 if ( !levelCombination1.compare("Multiplicative") || !levelCombination2.compare("Multiplicative") ) {
562
563 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Multiplicative Level Combination is not supported for NOX.");
564 ParameterListPtr_Type solverPList = this->getLinearSolverBuilder()->getNonconstParameterList();
565
566// pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",true);
567//
568// Teuchos::RCP<const Thyra::LinearOpBase<SC> > thyra_linOp = this->getPreconditioner()->getThyraPrec()->getUnspecifiedPrecOp();
569// Thyra::apply( *thyra_linOp, Thyra::NOTRANS, *thyraB, thyraX.ptr() );
570// pListThyraSolver->sublist("Preconditioner Types").sublist("FROSch").set("Only apply coarse",false);
571
572
573 }
574 }
575 }
576}
577#endif
578
579template<class SC,class LO,class GO,class NO>
580void NavierStokesAssFE<SC,LO,GO,NO>::calculateNonLinResidualVecWithMeshVelo(std::string type, double time, MultiVectorPtr_Type u_minus_w, MatrixPtr_Type P) const{
581
582
583 // this->reAssembleFSI( "FixedPoint", u_minus_w, P );
584
585 // We need to account for different parameters of time discretizations here
586 // This is ok for bdf with 1.0 scaling of the system. Would be wrong for Crank-Nicolson
587
588 this->system_->apply( *this->solution_, *this->residualVec_ );
589// this->residualVec_->getBlock(0)->writeMM("Ax.mm");
590// this->rhs_->getBlock(0)->writeMM("nsRHS.mm");
591 if (!type.compare("standard")){
592 this->residualVec_->update(-1.,*this->rhs_,1.);
593 if ( !this->sourceTerm_.is_null() )
594 this->residualVec_->update(-1.,*this->sourceTerm_,1.);
595 }
596 else if(!type.compare("reverse")){
597 this->residualVec_->update(1.,*this->rhs_,-1.); // this = -1*this + 1*rhs
598 if ( !this->sourceTerm_.is_null() )
599 this->residualVec_->update(1.,*this->sourceTerm_,1.);
600 }
601
602 // this might be set again by the TimeProblem after addition of M*u
603 this->bcFactory_->setBCMinusVector( this->residualVec_, this->solution_, time );
604
605// this->residualVec_->getBlock(0)->writeMM("b_Ax.mm");
606
607}
608
609
610template<class SC,class LO,class GO,class NO>
611Teuchos::RCP<Thyra::LinearOpBase<SC> > NavierStokesAssFE<SC,LO,GO,NO>::create_W_op() const
612{
613 //this->reAssemble("FixedPoint");
614 //this->reAssemble("Newton");
615
616 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
617 if ( !type.compare("Monolithic"))
618 return create_W_op_Monolithic( );
619 else if ( !type.compare("Teko")){
620#ifdef FEDD_HAVE_TEKO
621 return create_W_op_Block( );
622#else
623 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Teko not found! Build Trilinos with Teko.");
624#endif
625 }
626 else
627 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Unkown preconditioner/solver type.");
628}
629
630template<class SC,class LO,class GO,class NO>
631Teuchos::RCP<Thyra::LinearOpBase<SC> > NavierStokesAssFE<SC,LO,GO,NO>::create_W_op_Monolithic() const
632{
633 Teuchos::RCP<const Thyra::LinearOpBase<SC> > W_opConst = this->system_->getThyraLinOp();
634 Teuchos::RCP<Thyra::LinearOpBase<SC> > W_op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> >(W_opConst);
635 return W_op;
636}
637
638#ifdef FEDD_HAVE_TEKO
639template<class SC,class LO,class GO,class NO>
640Teuchos::RCP<Thyra::LinearOpBase<SC> > NavierStokesAssFE<SC,LO,GO,NO>::create_W_op_Block() const
641{
642
643 Teko::LinearOp thyraF = this->system_->getBlock(0,0)->getThyraLinOp();
644 Teko::LinearOp thyraBT = this->system_->getBlock(0,1)->getThyraLinOp();
645 Teko::LinearOp thyraB = this->system_->getBlock(1,0)->getThyraLinOp();
646
647 if (!this->system_->blockExists(1,1)){
648 MatrixPtr_Type dummy = Teuchos::rcp( new Matrix_Type( this->system_->getBlock(1,0)->getMap(), 1 ) );
649 dummy->fillComplete();
650 this->system_->addBlock( dummy, 1, 1 );
651 }
652
653 Teko::LinearOp thyraC = this->system_->getBlock(1,1)->getThyraLinOp();
654
655 Teuchos::RCP<const Thyra::LinearOpBase<SC> > W_opConst = Thyra::block2x2(thyraF,thyraBT,thyraB,thyraC);
656 Teuchos::RCP<Thyra::LinearOpBase<SC> > W_op = Teuchos::rcp_const_cast<Thyra::LinearOpBase<SC> >(W_opConst);
657 return W_op;
658}
659#endif
660
661template<class SC,class LO,class GO,class NO>
662Teuchos::RCP<Thyra::PreconditionerBase<SC> > NavierStokesAssFE<SC,LO,GO,NO>::create_W_prec() const
663{
664
665 this->initializeSolverBuilder();
666
667 std::string type = this->parameterList_->sublist("General").get("Preconditioner Method","Monolithic");
668 this->setBoundariesSystem();
669
670 if (!type.compare("Teko")) { //
671 this->setupPreconditioner( type );
672 stokesTekoPrecUsed_ = false;
673 }
674 else{
675 this->setupPreconditioner( type );
676 }
677
678 Teuchos::RCP<const Thyra::PreconditionerBase<SC> > thyraPrec = this->getPreconditionerConst()->getThyraPrecConst();
679 Teuchos::RCP<Thyra::PreconditionerBase<SC> > thyraPrecNonConst = Teuchos::rcp_const_cast<Thyra::PreconditionerBase<SC> >(thyraPrec);
680
681 return thyraPrecNonConst;
682
683}
684/*template<class SC,class LO,class GO,class NO>
685void NavierStokesAssFE<SC,LO,GO,NO>::reAssembleFSI(std::string type, MultiVectorPtr_Type u_minus_w, MatrixPtr_Type P) const {
686
687 if (this->verbose_)
688 std::cout << "-- Reassembly Navier-Stokes ("<< type <<") for FSI ... " << std::flush;
689
690 double density = this->parameterList_->sublist("Parameter").get("Density",1.);
691
692 MatrixPtr_Type ANW = Teuchos::rcp(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
693 if (type=="FixedPoint") {
694
695 MatrixPtr_Type N = Teuchos::rcp(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
696 this->feFactory_->assemblyAdvectionVecField( this->dim_, this->domain_FEType_vec_.at(0), N, u_minus_w, true );
697
698 N->resumeFill();
699 N->scale(density);
700 N->fillComplete( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getMapVecFieldUnique());
701 A_->addMatrix(1.,ANW,0.);
702
703 N->addMatrix(1.,ANW,1.);
704 // P must be scaled correctly in FSI
705 P->addMatrix(1.,ANW,1.);
706
707
708 }
709 else if(type=="Newton"){
710 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "reAssembleFSI should only be called for FPI-System.");
711 }
712 ANW->fillComplete( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getMapVecFieldUnique() );
713
714 this->system_->addBlock( ANW, 0, 0 );
715
716 if (this->verbose_)
717 std::cout << "done -- " << std::endl;
718}*/
719
720
721template<class SC,class LO,class GO,class NO>
722void NavierStokesAssFE<SC,LO,GO,NO>::reAssembleExtrapolation(BlockMultiVectorPtrArray_Type previousSolutions){
723
724 if (this->verbose_)
725 std::cout << "-- Reassembly Navier-Stokes (Extrapolation) ... " << std::flush;
726
727
728 double density = this->parameterList_->sublist("Parameter").get("Density",1.);
729
730 if (previousSolutions.size()>=2) {
731
732 MultiVectorPtr_Type extrapolatedVector = Teuchos::rcp( new MultiVector_Type( previousSolutions[0]->getBlock(0) ) );
733
734 extrapolatedVector->update( -1., *previousSolutions[1]->getBlock(0), 2. );
735
736 u_rep_->importFromVector(extrapolatedVector, true);
737 }
738 else if(previousSolutions.size()==1){
739 MultiVectorConstPtr_Type u = previousSolutions[0]->getBlock(0);
740 u_rep_->importFromVector(u, true);
741 }
742 else if (previousSolutions.size()==0){
743 MultiVectorConstPtr_Type u = this->solution_->getBlock(0);
744 u_rep_->importFromVector(u, true);
745 }
746
747 MatrixPtr_Type ANW = Teuchos::rcp(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
748
749 MatrixPtr_Type N = Teuchos::rcp(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
750 this->feFactory_->assemblyAdvectionVecField( this->dim_, this->domain_FEType_vec_.at(0), N, u_rep_, true );
751
752 N->resumeFill();
753 N->scale(density);
754 N->fillComplete( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getMapVecFieldUnique());
755
756 A_->addMatrix(1.,ANW,0.);
757 N->addMatrix(1.,ANW,1.);
758
759 ANW->fillComplete( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getMapVecFieldUnique() );
760
761 this->system_->addBlock( ANW, 0, 0 );
762
763 if (this->verbose_)
764 std::cout << "done -- " << std::endl;
765}
766
767
768
769
770
771
772// Use converged solution to compute viscosity field
773template<class SC,class LO,class GO,class NO>
774void NavierStokesAssFE<SC,LO,GO,NO>::computeSteadyPostprocessingViscosity_Solution() {
775
776
777 MultiVectorConstPtr_Type u = this->solution_->getBlock(0); // solution_ is initialized in problem_def.hpp so the most general class
778 u_rep_->importFromVector(u, true); // this is the current velocity solution at the nodes - distributed at the processors with repeated values
779
780 MultiVectorConstPtr_Type p = this->solution_->getBlock(1);
781 p_rep_->importFromVector(p, true); // this is the current pressure solution at the nodes - distributed at the processors with repeated values
782
783 // Reset here the viscosity so at this moment this makes only sense to call at the end of a simulation
784 // to visualize viscosity field
785 // @ToDo Add possibility for transient problem to save viscosity solution in each time step
786 viscosity_element_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getElementMap() ) );
787 this->feFactory_->computeSteadyViscosityFE_CM(this->dim_, this->getDomain(0)->getFEType(), this->getDomain(1)->getFEType(), this->dim_,1,u_rep_,p_rep_,this->parameterList_);
788
789 Teuchos::RCP<const MultiVector<SC,LO,GO,NO>> exportSolutionViscosityAssFE = this->feFactory_->const_output_fields->getBlock(0); // For now we assume that viscosity is always saved in first block
790 viscosity_element_->importFromVector(exportSolutionViscosityAssFE, true);
791
792
793
794}
795
796
797
798
799
800
801
802}
803
804#endif
FEDD::NavierStokesAssFE::calculateNonLinResidualVecWithMeshVelo
void calculateNonLinResidualVecWithMeshVelo(std::string type, double time, MultiVectorPtr_Type u_minus_w, MatrixPtr_Type P) const
Block Approach for Nonlinear Solver NOX. Input. Includes calculation of the residual vector and updat...
Definition NavierStokesAssFE_def.hpp:580
FEDD::NonLinearProblem
Definition NonLinearProblem_decl.hpp:24
FEDD::Problem< default_sc, default_lo, default_go, default_no >::sourceTerm_
BlockMultiVectorPtr_Type sourceTerm_
Definition Problem_decl.hpp:237
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement.cpp:5