Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
NonLinTPM_def.hpp
1#ifndef NonLinTPM_def_hpp
2#define NonLinTPM_def_hpp
3#include "NonLinTPM_decl.hpp"
4
13
14namespace FEDD {
15
16
17template<class SC,class LO,class GO,class NO>
18NonLinTPM<SC,LO,GO,NO>::NonLinTPM(const DomainConstPtr_Type &domainVelocity, std::string FETypeVelocity, const DomainConstPtr_Type &domainPressure, std::string FETypePressure, ParameterListPtr_Type parameterList):
19NonLinearProblem<SC,LO,GO,NO>(parameterList, domainVelocity->getComm())
20{
21
22 this->addVariable( domainVelocity , FETypeVelocity , "u" , domainVelocity->getDimension());
23 this->addVariable( domainPressure , FETypePressure , "p" , 1);
24 this->dim_ = this->getDomain(0)->getDimension();
25}
26
27template<class SC,class LO,class GO,class NO>
28NonLinTPM<SC,LO,GO,NO>::~NonLinTPM(){
29
30}
31
32template<class SC,class LO,class GO,class NO>
33void NonLinTPM<SC,LO,GO,NO>::info(){
34 this->infoProblem();
35}
36
37
38//template<class SC,class LO,class GO,class NO>
39//void NonLinTPM<SC,LO,GO,NO>::assemble( std::string type ){
40// if (type == "") {
41// if (this->verbose_)
42// std::cout << "-- Assembly ... " << std::flush;
43//
44// std::string tpmType = this->getParameterList()->sublist("Parameter").get("TPM Type","Biot");
45//
46// MapConstPtr_Type mapRepeatedConst1 = this->getDomain(0)->getMapRepeated();
47// MapConstPtr_Type mapRepeatedConst2 = this->getDomain(1)->getMapRepeated();
48// MapPtr_Type mapRepeated1 = Teuchos::rcp_const_cast<Map_Type>(mapRepeatedConst1);
49// MapPtr_Type mapRepeated2 = Teuchos::rcp_const_cast<Map_Type>(mapRepeatedConst2);
50//
51// MatrixPtr_Type A00(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
52// MatrixPtr_Type A01(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(1)->getDimension() * this->getDomain(1)->getApproxEntriesPerRow() ) );
53// MatrixPtr_Type A10(new Matrix_Type( this->getDomain(1)->getMapUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
54// MatrixPtr_Type A11(new Matrix_Type( this->getDomain(1)->getMapUnique(), this->getDomain(1)->getDimension() * this->getDomain(1)->getApproxEntriesPerRow() ) );
55//
56// MultiVectorPtr_Type F0(new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated(), 1 ) );
57// MultiVectorPtr_Type F1(new MultiVector_Type( this->getDomain(1)->getMapRepeated(), 1 ) );
58//
59// if (u_repNewton_.is_null()){
60// u_repNewton_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated() )) ;
61// u_repNewton_->putScalar(0.);
62// }
63// if (p_repNewton_.is_null()){
64// p_repNewton_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(1)->getMapRepeated() )) ;
65// p_repNewton_->putScalar(0.);
66// }
67//
68// if (u_repTime_.is_null()){
69// u_repTime_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated() )) ;
70// u_repTime_->putScalar(0.);
71// }
72// if (p_repTime_.is_null()){
73// p_repTime_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(1)->getMapRepeated() )) ;
74// p_repTime_->putScalar(0.);
75// }
76// if(tpmType == "Biot" || tpmType == "Biot-StVK"){
77// this->feFactory_->assemblyAceGenTPM(A00,
78// A01,
79// A10,
80// A11,
81// F0,
82// F1,
83// mapRepeated1,
84// mapRepeated2,
85// this->parameterList_,
86// u_repNewton_,
87// p_repNewton_,
88// u_repTime_,
89// p_repTime_,
90// true,
91// false);
92// }
93//
94// this->system_.reset(new BlockMatrix_Type(2));
95// this->system_->addBlock( A00, 0, 0 );
96// this->system_->addBlock( A01, 0, 1 );
97// this->system_->addBlock( A10, 1, 0 );
98// this->system_->addBlock( A11, 1, 1 );
99//
100// this->initializeVectors();
101// this->initializeVectorsNonLinear();
102//
103// MultiVectorPtr_Type F0Unique = Teuchos::rcp(new MultiVector_Type( this->getDomain(0)->getMapVecFieldUnique() ) );
104// MultiVectorPtr_Type F1Unique = Teuchos::rcp(new MultiVector_Type( this->getDomain(1)->getMapUnique() ) );
105// F0Unique->exportFromVector( F0, false, "Add" );
106// F1Unique->exportFromVector( F1, false, "Add" );
107//
108// this->rhs_->addBlock( F0Unique, 0 );
109// this->rhs_->addBlock( F1Unique, 1 );
110//
111// this->assembleSourceTerm( 0. );
112// this->addToRhs( this->sourceTerm_ );
113//
114// if (this->verbose_)
115// std::cout << "done -- " << std::endl;
116// } else {
117// this->reAssmble(type);
118// }
119//}
120
121
122template<class SC,class LO,class GO,class NO>
123void NonLinTPM<SC,LO,GO,NO>::calculateNonLinResidualVec(std::string type, double time) const{
124
125 this->assemble("Newton-Residual");
126 if (type == "external" ){// we assume reverse computation: Ax-b but we need to scale the sourceTerm with -1
127 this->residualVec_->update(1.,*this->rhs_,0.);
128 if ( !this->sourceTerm_.is_null() ){
129 this->residualVec_->update(-1.,*this->sourceTerm_,1.);
130 }
131 }
132 else{
133 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "Unknown type for residual computation.");
134 }
135 this->bcFactory_->setBCMinusVector( this->residualVec_, this->solution_, time );
136
137}
138
139template<class SC,class LO,class GO,class NO>
140void NonLinTPM<SC,LO,GO,NO>::reAssemble( BlockMultiVectorPtr_Type previousSolution ) const {
141 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error, "This should not be used.");
142 MultiVectorConstPtr_Type u_prev = previousSolution->getBlock(0);
143 MultiVectorConstPtr_Type p_prev = previousSolution->getBlock(1);
144 u_repTime_->importFromVector(u_prev, true);
145 p_repTime_->importFromVector(p_prev, true);
146}
147
148
149template<class SC,class LO,class GO,class NO>
150void NonLinTPM<SC,LO,GO,NO>::assemble( std::string type ) const{
151
152 if (this->verbose_)
153 std::cout << "-- Assembly nonlinear TPM " <<" ("<<type <<") ... " << std::flush;
154
155 std::string tpmType = this->parameterList_->sublist("Parameter").get("TPM Type","Biot");
156
157 if (type=="Newton-Residual") {
158 type = "Assemble";
159 }
160
161 if (type == "FirstAssemble") {
162 if (u_repNewton_.is_null()){
163 u_repNewton_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated() )) ;
164 u_repNewton_->putScalar(0.);
165 }
166 if (p_repNewton_.is_null()){
167 p_repNewton_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(1)->getMapRepeated() )) ;
168 p_repNewton_->putScalar(0.);
169 }
170
171 if (u_repTime_.is_null()){
172 u_repTime_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated() )) ;
173 u_repTime_->putScalar(0.);
174 }
175 if (p_repTime_.is_null()){
176 p_repTime_ = Teuchos::rcp( new MultiVector_Type( this->getDomain(1)->getMapRepeated() )) ;
177 p_repTime_->putScalar(0.);
178 }
179 }
180
181
182 if (type == "SetSolutionInTime") {
183 TEUCHOS_TEST_FOR_EXCEPTION( u_repTime_.is_null(), std::runtime_error, "u_repTime_ not initialized.");
184 TEUCHOS_TEST_FOR_EXCEPTION( p_repTime_.is_null(), std::runtime_error, "p_repTime_ not initialized.");
185 MultiVectorConstPtr_Type u = this->solution_->getBlock(0);
186 MultiVectorConstPtr_Type p = this->solution_->getBlock(1);
187 u_repTime_->importFromVector(u, true);
188 p_repTime_->importFromVector(p, true);
189 }
190 else if (type == "SetSolutionNewton") {
191 TEUCHOS_TEST_FOR_EXCEPTION( u_repNewton_.is_null(), std::runtime_error, "u_repNewton_ not initialized.");
192 TEUCHOS_TEST_FOR_EXCEPTION( p_repNewton_.is_null(), std::runtime_error, "p_repNewton_ not initialized.");
193
194 MultiVectorConstPtr_Type u = this->solution_->getBlock(0);
195 MultiVectorConstPtr_Type p = this->solution_->getBlock(1);
196 u_repNewton_->importFromVector(u, true);
197 p_repNewton_->importFromVector(p, true);
198 }
199 else if( type == "FirstAssemble" || type == "Assemble" || type == "OnlyUpdate" || type == "AssembleAndUpdate"){
200 if (this->verbose_)
201 std::cout << "-- External assembly ... " << std::flush;
202
203 MapConstPtr_Type mapRepeatedConst1 = this->getDomain(0)->getMapRepeated();
204 MapConstPtr_Type mapRepeatedConst2 = this->getDomain(1)->getMapRepeated();
205 MapPtr_Type mapRepeated1 = Teuchos::rcp_const_cast<Map_Type>(mapRepeatedConst1);
206 MapPtr_Type mapRepeated2 = Teuchos::rcp_const_cast<Map_Type>(mapRepeatedConst2);
207
208 MatrixPtr_Type A00(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
209 MatrixPtr_Type A01(new Matrix_Type( this->getDomain(0)->getMapVecFieldUnique(), this->getDomain(1)->getDimension() * this->getDomain(1)->getApproxEntriesPerRow() ) );
210 MatrixPtr_Type A10(new Matrix_Type( this->getDomain(1)->getMapUnique(), this->getDomain(0)->getDimension() * this->getDomain(0)->getApproxEntriesPerRow() ) );
211 MatrixPtr_Type A11(new Matrix_Type( this->getDomain(1)->getMapUnique(), this->getDomain(1)->getDimension() * this->getDomain(1)->getApproxEntriesPerRow() ) );
212
213
214 MultiVectorPtr_Type F0(new MultiVector_Type( this->getDomain(0)->getMapVecFieldRepeated(), 1 ) );
215 MultiVectorPtr_Type F1(new MultiVector_Type( this->getDomain(1)->getMapRepeated(), 1 ) );
216
217 bool update(type == "FirstAssemble" || type == "Assemble" || type == "AssembleAndUpdate");
218 bool updateHistory(type == "OnlyUpdate" || type == "AssembleAndUpdate");
219
220 if(tpmType == "Biot" || tpmType == "Biot-StVK"){
221 this->feFactory_->assemblyAceGenTPM(A00,
222 A01,
223 A10,
224 A11,
225 F0,
226 F1,
227 mapRepeated1,
228 mapRepeated2,
229 this->parameterList_,
230 u_repNewton_,
231 p_repNewton_,
232 u_repTime_,
233 p_repTime_,
234 update,
235 updateHistory);
236 }
237 if (update) {
238 this->system_->addBlock( A00, 0, 0 );
239 this->system_->addBlock( A01, 0, 1 );
240 this->system_->addBlock( A10, 1, 0 );
241 this->system_->addBlock( A11, 1, 1 );
242
243
244 MultiVectorPtr_Type F0Unique = Teuchos::rcp(new MultiVector_Type( this->getDomain(0)->getMapVecFieldUnique() ) );
245 MultiVectorPtr_Type F1Unique = Teuchos::rcp(new MultiVector_Type( this->getDomain(1)->getMapUnique() ) );
246 F0Unique->exportFromVector( F0, false, "Add" );
247 F1Unique->exportFromVector( F1, false, "Add" );
248
249 this->rhs_->addBlock( F0Unique, 0 );
250 this->rhs_->addBlock( F1Unique, 1 );
251 }
252 }
253 if (this->verbose_)
254 std::cout << "done -- " << std::endl;
255
256}
257
258template<class SC,class LO,class GO,class NO>
259void NonLinTPM<SC,LO,GO,NO>::reAssemble( MatrixPtr_Type& massmatrix, std::string type ) const
260{
261
262}
263
264template<class SC,class LO,class GO,class NO>
265void NonLinTPM<SC,LO,GO,NO>::reAssembleExtrapolation(BlockMultiVectorPtrArray_Type previousSolutions){
266
267
268 TEUCHOS_TEST_FOR_EXCEPTION( true, std::logic_error, "Only Newton for NonLinTPM!");
269
270}
271
272template<class SC,class LO,class GO,class NO>
273void NonLinTPM<SC,LO,GO,NO>::evalModelImpl(const Thyra::ModelEvaluatorBase::InArgs<SC> &inArgs,
274 const Thyra::ModelEvaluatorBase::OutArgs<SC> &outArgs
275 ) const
276{
277 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented! NonLinTPM evalModelImpl(...)");
278
279}
280
281template<class SC,class LO,class GO,class NO>
282Teuchos::RCP<Thyra::LinearOpBase<SC> > NonLinTPM<SC,LO,GO,NO>::create_W_op() const
283{
284 Teuchos::RCP<Thyra::LinearOpBase<SC> > dummy;
285 return dummy;
286}
287
288template<class SC,class LO,class GO,class NO>
289Teuchos::RCP<Thyra::PreconditionerBase<SC> > NonLinTPM<SC,LO,GO,NO>::create_W_prec() const
290{
291 Teuchos::RCP<Thyra::PreconditionerBase<SC> > dummy;
292 return dummy;
293}
294
295}
296#endif
FEDD::NonLinTPM::assemble
virtual void assemble(std::string type="") const
assemble of type exectuted by the derived specific non-linear problem classes
Definition NonLinTPM_def.hpp:150
FEDD::NonLinTPM::calculateNonLinResidualVec
virtual void calculateNonLinResidualVec(std::string type, double time=0.) const
Virtual function which is implemented in the specific non-linear problem classes to calculate the non...
Definition NonLinTPM_def.hpp:123
FEDD::NonLinearProblem
Definition NonLinearProblem_decl.hpp:28
FEDD::Problem< default_sc, default_lo, default_go, default_no >::sourceTerm_
BlockMultiVectorPtr_Type sourceTerm_
Definition Problem_decl.hpp:251
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36