Panzer Version of the Day
Loading...
Searching...
No Matches
Panzer_InitialCondition_Builder.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Panzer: A partial differential equation assembly
5// engine for strongly coupled complex multiphysics systems
6// Copyright (2011) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roger P. Pawlowski (rppawlo@sandia.gov) and
39// Eric C. Cyr (eccyr@sandia.gov)
40// ***********************************************************************
41// @HEADER
42
44
45#include "Teuchos_Assert.hpp"
46
47#include "Panzer_EquationSet_DefaultImpl.hpp"
52
53namespace panzer {
54
55void
57 const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
59 const Teuchos::ParameterList& ic_block_closure_models,
61 const Teuchos::ParameterList& user_data,
62 const bool write_graphviz_file,
63 const std::string& graphviz_file_prefix,
64 std::map< std::string, Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >& phx_ic_field_managers)
65{
66 std::vector<Teuchos::RCP<panzer::PhysicsBlock> >::const_iterator blkItr;
67 for (blkItr=physicsBlocks.begin();blkItr!=physicsBlocks.end();++blkItr) {
68 Teuchos::RCP<panzer::PhysicsBlock> pb = *blkItr;
69 std::string blockId = pb->elementBlockID();
70
71 // build a field manager object
72 Teuchos::RCP<PHX::FieldManager<panzer::Traits> > fm
73 = Teuchos::rcp(new PHX::FieldManager<panzer::Traits>);
74
75 // Choose model sublist for this element block
76 std::string closure_model_name = "";
77 if (ic_block_closure_models.isSublist(blockId))
78 closure_model_name = blockId;
79 else if (ic_block_closure_models.isSublist("Default"))
80 closure_model_name = "Default";
81 else
82 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Failed to find initial condition for element block \"" << blockId
83 << "\". You must provide an initial condition for each element block or set a default!"
84 << ic_block_closure_models);
85
86 // Only the residual is used by initial conditions
87 std::vector<bool> active_evaluation_types(Sacado::mpl::size<panzer::Traits::EvalTypes>::value,false);
88 int residual_index = Sacado::mpl::find<panzer::Traits::EvalTypes,panzer::Traits::Residual>::value;
89 active_evaluation_types[residual_index] = true;
90 pb->setActiveEvaluationTypes(active_evaluation_types);
91
92 // use the physics block to register evaluators
93 pb->buildAndRegisterInitialConditionEvaluators(*fm, cm_factory, closure_model_name, ic_block_closure_models, lo_factory, user_data);
94
95 pb->activateAllEvaluationTypes();
96
97 // build the setup data using passed in information
98 Traits::SD setupData;
99 const WorksetDescriptor wd = blockDescriptor(blockId);
100 setupData.worksets_ = wkstContainer.getWorksets(wd);
101 setupData.orientations_ = wkstContainer.getOrientations();
102
103 int i=0;
104 for (auto eval_type=fm->begin(); eval_type != fm->end(); ++eval_type,++i) {
105 if (active_evaluation_types[i])
106 eval_type->postRegistrationSetup(setupData,*fm,false,false,nullptr);
107 }
108
109 phx_ic_field_managers[blockId] = fm;
110
111 if (write_graphviz_file)
112 fm->writeGraphvizFile(graphviz_file_prefix+"_IC_"+blockId);
113 }
114}
115
116void
118 const std::vector<Teuchos::RCP<panzer::PhysicsBlock> >& physicsBlocks,
120 const Teuchos::ParameterList& closure_models,
121 const Teuchos::ParameterList& ic_block_closure_models,
123 const Teuchos::ParameterList& user_data,
124 const bool write_graphviz_file,
125 const std::string& graphviz_file_prefix,
126 std::map< std::string, Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >& phx_ic_field_managers)
127{
128 std::vector<Teuchos::RCP<panzer::PhysicsBlock> >::const_iterator blkItr;
129 for (blkItr=physicsBlocks.begin();blkItr!=physicsBlocks.end();++blkItr) {
130 Teuchos::RCP<panzer::PhysicsBlock> pb = *blkItr;
131 std::string blockId = pb->elementBlockID();
132
133 // build a field manager object
134 Teuchos::RCP<PHX::FieldManager<panzer::Traits> > fm
135 = Teuchos::rcp(new PHX::FieldManager<panzer::Traits>);
136
137 // Choose model sublist for this element block
138 std::string closure_model_name = "";
139 if (ic_block_closure_models.isSublist(blockId))
140 closure_model_name = blockId;
141 else if (ic_block_closure_models.isSublist("Default"))
142 closure_model_name = "Default";
143 else
144 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Failed to find initial condition for element block \"" << blockId
145 << "\". You must provide an initial condition for each element block or set a default!"
146 << ic_block_closure_models);
147
148 // Only the residual is used by initial conditions
149 std::vector<bool> active_evaluation_types(Sacado::mpl::size<panzer::Traits::EvalTypes>::value,false);
150 int residual_index = Sacado::mpl::find<panzer::Traits::EvalTypes,panzer::Traits::Residual>::value;
151 active_evaluation_types[residual_index] = true;
152 pb->setActiveEvaluationTypes(active_evaluation_types);
153
154 // build and register all closure models
155 pb->buildAndRegisterClosureModelEvaluators(*fm,cm_factory,closure_models,user_data);
156
157 // use the physics block to register evaluators
158 pb->buildAndRegisterInitialConditionEvaluators(*fm, cm_factory, closure_model_name, ic_block_closure_models, lo_factory, user_data);
159
160 pb->activateAllEvaluationTypes();
161
162 // build the setup data using passed in information
163 Traits::SD setupData;
164 const WorksetDescriptor wd = blockDescriptor(blockId);
165 setupData.worksets_ = wkstContainer.getWorksets(wd);
166 setupData.orientations_ = wkstContainer.getOrientations();
167
168 int i=0;
169 for (auto eval_type=fm->begin(); eval_type != fm->end(); ++eval_type,++i) {
170 if (active_evaluation_types[i])
171 eval_type->postRegistrationSetup(setupData,*fm,false,false,nullptr);
172 }
173
174 phx_ic_field_managers[blockId] = fm;
175
176 if (write_graphviz_file)
177 fm->writeGraphvizFile(graphviz_file_prefix+"_IC_"+blockId);
178 }
179}
180
181void
183 const std::map< std::string,Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >& phx_ic_field_managers,
184 Teuchos::RCP<panzer::LinearObjContainer> loc,
186 const double time_stamp)
187{
188 typedef LinearObjContainer LOC;
190
191 // allocate a ghosted container for the initial condition
192 Teuchos::RCP<LOC> ghostedloc = lo_factory.buildGhostedLinearObjContainer();
193 lo_factory.initializeGhostedContainer(LOC::X,*ghostedloc);
194
195 // allocate a counter to keep track of where this processor set initial conditions
196 Teuchos::RCP<LOC> localCounter = lo_factory.buildPrimitiveGhostedLinearObjContainer();
197 Teuchos::RCP<LOC> globalCounter = lo_factory.buildPrimitiveLinearObjContainer();
198 Teuchos::RCP<LOC> summedGhostedCounter = lo_factory.buildPrimitiveGhostedLinearObjContainer();
199
200 lo_factory.initializeGhostedContainer(LOC::F,*localCounter); // store counter in F
201 localCounter->initialize();
202
203 ped.gedc->addDataObject("Residual Scatter Container",ghostedloc);
204 ped.gedc->addDataObject("Dirichlet Counter",localCounter);
206
207 for(std::map< std::string,Teuchos::RCP< PHX::FieldManager<panzer::Traits> > >::const_iterator itr=phx_ic_field_managers.begin();
208 itr!=phx_ic_field_managers.end();++itr) {
209 std::string blockId = itr->first;
210 Teuchos::RCP< PHX::FieldManager<panzer::Traits> > fm = itr->second;
211
212 fm->preEvaluate<panzer::Traits::Residual>(ped);
213
214 // Loop over worksets in this element block
215 const WorksetDescriptor wd = blockDescriptor(blockId);
216 std::vector<panzer::Workset>& w = *wkstContainer.getWorksets(wd);
217 for (std::size_t i = 0; i < w.size(); ++i) {
218 panzer::Workset& workset = w[i];
219
220 // Need to figure out how to get restart time from Rythmos.
221 workset.time = time_stamp;
222
223 fm->evaluateFields<panzer::Traits::Residual>(workset);
224 }
225 }
226
227 lo_factory.initializeGhostedContainer(LOC::F,*summedGhostedCounter); // store counter in F
228 summedGhostedCounter->initialize();
229
230 // do communication to build summed ghosted counter for dirichlet conditions
231 {
232 lo_factory.initializeContainer(LOC::F,*globalCounter); // store counter in F
233 globalCounter->initialize();
234 lo_factory.ghostToGlobalContainer(*localCounter,*globalCounter,LOC::F);
235 // Here we do the reduction across all processors so that the number of times
236 // a dirichlet condition is applied is summed into the global counter
237
238 lo_factory.globalToGhostContainer(*globalCounter,*summedGhostedCounter,LOC::F);
239 // finally we move the summed global vector into a local ghosted vector
240 // so that the dirichlet conditions can be applied to both the ghosted
241 // right hand side and the ghosted matrix
242 }
243
245 gedc.addDataObject("Residual Scatter Container",ghostedloc);
246
247 // adjust ghosted system for boundary conditions
248 for(GlobalEvaluationDataContainer::iterator itr=gedc.begin();itr!=gedc.end();itr++) {
249 Teuchos::RCP<LOC> loc2 = Teuchos::rcp_dynamic_cast<LOC>(itr->second);
250 if(loc2!=Teuchos::null) {
251 bool zeroVectorRows = false;
252 bool adjustX = true;
253 lo_factory.adjustForDirichletConditions(*localCounter,*summedGhostedCounter,*loc2, zeroVectorRows, adjustX);
254 }
255 }
256
257 // gather the ghosted solution back into the global container, which performs the sum
258 lo_factory.ghostToGlobalContainer(*ghostedloc,*loc,LOC::X);
259}
260
261// This is an anonymous namespace that hides a few helper classes for the control
262// initial condition construction. In particual an IC Equation set is defined that
263// is useful for building initial condition vectors.
264namespace {
265
266template <typename EvalT>
267class EquationSet_IC : public EquationSet_DefaultImpl<EvalT> {
268public:
269
273 EquationSet_IC(const Teuchos::RCP<Teuchos::ParameterList>& params,
274 const int& default_integration_order,
275 const CellData& cell_data,
276 const Teuchos::RCP<GlobalData>& global_data,
277 const bool build_transient_support);
278
281 void buildAndRegisterEquationSetEvaluators(PHX::FieldManager<Traits>& /* fm */,
282 const FieldLibrary& /* field_library */,
283 const Teuchos::ParameterList& /* user_data */) const {}
284
285};
286
287// ***********************************************************************
288template <typename EvalT>
289EquationSet_IC<EvalT>::
290EquationSet_IC(const Teuchos::RCP<Teuchos::ParameterList>& params,
291 const int& default_integration_order,
292 const CellData& cell_data,
293 const Teuchos::RCP<GlobalData>& global_data,
294 const bool build_transient_support) :
295 EquationSet_DefaultImpl<EvalT>(params, default_integration_order, cell_data, global_data, build_transient_support )
296{
297 // ********************
298 // Validate and parse parameter list
299 // ********************
300 Teuchos::ParameterList valid_parameters_sublist;
301 valid_parameters_sublist.set("Basis Type","HGrad","Type of Basis to use");
302 valid_parameters_sublist.set("Basis Order",1,"Order of the basis");
303
304 for(auto itr=params->begin();itr!=params->end();++itr) {
305
306 const std::string field = params->name(itr);
307 const Teuchos::ParameterEntry & entry = params->entry(itr);
308
309 // only process lists
310 if(!entry.isList()) continue;
311
312 Teuchos::ParameterList & basisPL = entry.getValue((Teuchos::ParameterList *) 0);
313 basisPL.validateParametersAndSetDefaults(valid_parameters_sublist);
314
315 std::string basis_type = basisPL.get<std::string>("Basis Type");
316 int basis_order = basisPL.get<int>("Basis Order");
317
318 this->addDOF(field,basis_type,basis_order,default_integration_order);
319 }
320
321 this->addClosureModel("");
322
323 this->setupDOFs();
324}
325
326PANZER_DECLARE_EQSET_TEMPLATE_BUILDER(EquationSet_IC, EquationSet_IC)
327
328// A user written factory that creates each equation set. The key member here
329// is buildEquationSet
330class IC_EquationSetFactory : public EquationSetFactory {
331public:
332
333 Teuchos::RCP<EquationSet_TemplateManager<Traits> >
334 buildEquationSet(const Teuchos::RCP<Teuchos::ParameterList>& params,
335 const int& default_integration_order,
336 const CellData& cell_data,
337 const Teuchos::RCP<GlobalData>& global_data,
338 const bool build_transient_support) const
339 {
340 Teuchos::RCP<EquationSet_TemplateManager<Traits> > eq_set=
341 Teuchos::rcp(new EquationSet_TemplateManager<Traits>);
342
343 bool found = false; // this is used by PANZER_BUILD_EQSET_OBJECTS
344
345 // macro checks if(ies.name=="Poisson") then an EquationSet_Energy object is constructed
346 PANZER_BUILD_EQSET_OBJECTS("IC", EquationSet_IC)
347
348 // make sure your equation set has been found
349 if(!found) {
350 std::string msg = "Error - the \"Equation Set\" called \"" + params->get<std::string>("Type") +
351 "\" is not a valid equation set identifier. Please supply the correct factory.\n";
352 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, msg);
353 }
354
355 return eq_set;
356 }
357
358};
359
360} // end anonymous namespace
361
362void
363setupControlInitialCondition(const std::map<std::string,Teuchos::RCP<const shards::CellTopology> > & block_ids_to_cell_topo,
364 const std::map<std::string,std::vector<ICFieldDescriptor> > & block_ids_to_fields,
365 WorksetContainer & wkstContainer,
366 const LinearObjFactory<Traits> & lof,
368 const Teuchos::ParameterList & ic_closure_models,
369 const Teuchos::ParameterList & user_data,
370 int workset_size,
371 double t0,
372 const Teuchos::RCP<Thyra::VectorBase<double> > & vec)
373{
374 std::vector<Teuchos::RCP<PhysicsBlock> > physics_blocks;
375 buildICPhysicsBlocks(block_ids_to_cell_topo,block_ids_to_fields,workset_size,physics_blocks);
376
377 std::map<std::string, Teuchos::RCP< PHX::FieldManager<Traits> > > phx_ic_field_managers;
379 physics_blocks,
380 cm_factory,
381 ic_closure_models,
382 lof,
383 user_data,
384 false,
385 "initial_condition_control_test",
386 phx_ic_field_managers);
387
388
389 Teuchos::RCP<LinearObjContainer> loc = lof.buildLinearObjContainer();
390 Teuchos::rcp_dynamic_cast<ThyraObjContainer<double> >(loc)->set_x_th(vec);
391
392 evaluateInitialCondition(wkstContainer, phx_ic_field_managers, loc, lof, t0);
393}
394
395
396void
397buildICPhysicsBlocks(const std::map<std::string,Teuchos::RCP<const shards::CellTopology> > & block_ids_to_cell_topo,
398 const std::map<std::string,std::vector<ICFieldDescriptor> > & block_ids_to_fields,
399 int workset_size,
400 std::vector<Teuchos::RCP<PhysicsBlock> > & physics_blocks)
401{
402 using Teuchos::RCP;
403 using Teuchos::rcp;
404
405 std::map<std::string,std::string> block_ids_to_physics_ids;
406
407 RCP<Teuchos::ParameterList> ipb = rcp(new Teuchos::ParameterList);
408
409 for(auto itr=block_ids_to_cell_topo.begin();itr!=block_ids_to_cell_topo.end();++itr) {
410 std::string eblock = itr->first;
411 RCP<const shards::CellTopology> ct = itr->second;
412
413 // get the field descriptor vector, check to make sure block is represented
414 auto fds_itr = block_ids_to_fields.find(eblock);
415 TEUCHOS_ASSERT(fds_itr!=block_ids_to_fields.end());
416
417 const std::vector<ICFieldDescriptor> & fd_vec = fds_itr->second;
418
419 std::string physics_id = "ic_"+eblock;
420 block_ids_to_physics_ids[eblock] = physics_id;
421
422 // start building a physics block named "ic_"+eblock, with one anonymous list
423 Teuchos::ParameterList & physics_block = ipb->sublist(physics_id).sublist("");
424 physics_block.set("Type","IC"); // point IC type
425
426 for(std::size_t i=0;i<fd_vec.size();i++) {
427 const ICFieldDescriptor & fd = fd_vec[i];
428
429 // number the lists, these should be anonymous
430 physics_block.sublist(fd.fieldName).set("Basis Type", fd.basisName);
431 physics_block.sublist(fd.fieldName).set("Basis Order",fd.basisOrder);
432 }
433 }
434
435 RCP<EquationSetFactory> eqset_factory = Teuchos::rcp(new IC_EquationSetFactory);
436
437 RCP<GlobalData> gd = createGlobalData();
438 buildPhysicsBlocks(block_ids_to_physics_ids,
439 block_ids_to_cell_topo,
440 ipb,1,workset_size,eqset_factory,gd,false,physics_blocks);
441}
442
443} // end namespace panzer
#define PANZER_BUILD_EQSET_OBJECTS(key, fType)
#define PANZER_DECLARE_EQSET_TEMPLATE_BUILDER(fClass, fType)
PHX::MDField< ScalarT, panzer::Cell, panzer::BASIS > field
A field to which we'll contribute, or in which we'll store, the result of computing this integral.
void addDataObject(const std::string &key, const Teuchos::RCP< GlobalEvaluationData > &ged)
std::unordered_map< std::string, Teuchos::RCP< GlobalEvaluationData > >::iterator iterator
virtual Teuchos::RCP< LinearObjContainer > buildLinearObjContainer() const =0
virtual void initializeGhostedContainer(int, LinearObjContainer &loc) const =0
virtual Teuchos::RCP< LinearObjContainer > buildGhostedLinearObjContainer() const =0
virtual Teuchos::RCP< LinearObjContainer > buildPrimitiveGhostedLinearObjContainer() const =0
virtual void globalToGhostContainer(const LinearObjContainer &container, LinearObjContainer &ghostContainer, int) const =0
virtual void ghostToGlobalContainer(const LinearObjContainer &ghostContainer, LinearObjContainer &container, int) const =0
virtual Teuchos::RCP< LinearObjContainer > buildPrimitiveLinearObjContainer() const =0
virtual void initializeContainer(int, LinearObjContainer &loc) const =0
virtual void adjustForDirichletConditions(const LinearObjContainer &localBCRows, const LinearObjContainer &globalBCRows, LinearObjContainer &ghostedObjs, bool zeroVectorRows=false, bool adjustX=false) const =0
Class that provides access to worksets on each element block and side set.
Teuchos::RCP< std::vector< Workset > > getWorksets(const WorksetDescriptor &wd)
Access to volume worksets.
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > getOrientations() const
WorksetDescriptor blockDescriptor(const std::string &eBlock)
void evaluateInitialCondition(WorksetContainer &wkstContainer, const std::map< std::string, Teuchos::RCP< PHX::FieldManager< panzer::Traits > > > &phx_ic_field_managers, Teuchos::RCP< panzer::LinearObjContainer > loc, const panzer::LinearObjFactory< panzer::Traits > &lo_factory, const double time_stamp)
void buildICPhysicsBlocks(const std::map< std::string, Teuchos::RCP< const shards::CellTopology > > &block_ids_to_cell_topo, const std::map< std::string, std::vector< ICFieldDescriptor > > &block_ids_to_fields, int workset_size, std::vector< Teuchos::RCP< PhysicsBlock > > &physics_blocks)
Teuchos::RCP< panzer::GlobalData > createGlobalData(bool build_default_os, int print_process)
void setupInitialConditionFieldManagers(WorksetContainer &wkstContainer, const std::vector< Teuchos::RCP< panzer::PhysicsBlock > > &physicsBlocks, const panzer::ClosureModelFactory_TemplateManager< panzer::Traits > &cm_factory, const Teuchos::ParameterList &ic_block_closure_models, const panzer::LinearObjFactory< panzer::Traits > &lo_factory, const Teuchos::ParameterList &user_data, const bool write_graphviz_file, const std::string &graphviz_file_prefix, std::map< std::string, Teuchos::RCP< PHX::FieldManager< panzer::Traits > > > &phx_ic_field_managers)
Builds PHX::FieldManager objects for inital conditions and registers evaluators.
void setupControlInitialCondition(const std::map< std::string, Teuchos::RCP< const shards::CellTopology > > &block_ids_to_cell_topo, const std::map< std::string, std::vector< ICFieldDescriptor > > &block_ids_to_fields, WorksetContainer &wkstContainer, const LinearObjFactory< Traits > &lof, const ClosureModelFactory_TemplateManager< Traits > &cm_factory, const Teuchos::ParameterList &ic_closure_models, const Teuchos::ParameterList &user_data, int workset_size, double t0, const Teuchos::RCP< Thyra::VectorBase< double > > &vec)
Teuchos::RCP< GlobalEvaluationDataContainer > gedc
std::string first_sensitivities_name
Teuchos::RCP< const std::vector< panzer::Workset > > worksets_
Teuchos::RCP< const std::vector< Intrepid2::Orientation > > orientations_