31 using namespace MoFEM;
36 #include <boost/numeric/ublas/vector_proxy.hpp>
37 #include <boost/numeric/ublas/matrix.hpp>
38 #include <boost/numeric/ublas/matrix_proxy.hpp>
39 #include <boost/numeric/ublas/vector.hpp>
40 #include <adolc/adolc.h>
45 #include <boost/program_options.hpp>
47 namespace po = boost::program_options;
49 #include <MethodForForceScaling.hpp>
60 #include <TimeForceScale.hpp>
62 #include <boost/iostreams/tee.hpp>
63 #include <boost/iostreams/stream.hpp>
66 namespace bio = boost::iostreams;
67 using bio::tee_device;
72 static char help[] =
"...\n\n";
81 int main(
int argc,
char *argv[]) {
90 PetscBool flg_gel_config,flg_file;
92 char gel_config_file[255];
94 PetscBool is_partitioned = PETSC_FALSE;
95 PetscBool is_atom_test = PETSC_FALSE;
99 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
"",
"Elastic Config",
"none"); CHKERRQ(
ierr);
100 ierr = PetscOptionsString(
108 ierr = PetscOptionsInt(
110 "default approximation order",
117 ierr = PetscOptionsBool(
118 "-my_is_partitioned",
119 "set if mesh is partitioned (this result that each process keeps only part of the mes",
"",
120 PETSC_FALSE,&is_partitioned,PETSC_NULL
122 ierr = PetscOptionsString(
124 "gel configuration file name",
"",
125 "gel_config.in",gel_config_file,255,&flg_gel_config
128 ierr = PetscOptionsBool(
130 "is used with testing, exit with error when diverged",
"",
131 PETSC_FALSE,&is_atom_test,PETSC_NULL
134 ierr = PetscOptionsEnd(); CHKERRQ(
ierr);
137 if(flg_file != PETSC_TRUE) {
138 SETERRQ(PETSC_COMM_SELF,1,
"*** ERROR -my_file (MESH FILE NEEDED)");
141 if(is_partitioned == PETSC_TRUE) {
145 "PARALLEL=BCAST_DELETE;"
146 "PARALLEL_RESOLVE_SHARED_ENTS;"
147 "PARTITION=PARALLEL_PARTITION;";
154 ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,
MYPCOMM_INDEX);
155 if(pcomm == NULL) pcomm =
new ParallelComm(&moab,PETSC_COMM_WORLD);
169 map<int,ThermalElement::FluxData> set_of_solvent_fluxes;
176 bool check_if_spatial_field_exist = m_field.
check_field(
"SPATIAL_POSITION");
236 map<int,BlockData> block_data;
239 ifstream ini_file(gel_config_file);
240 po::variables_map vm;
241 po::options_description config_file_options;
245 rval = moab.tag_get_handle(
246 "BLOCK_ID",1,MB_TYPE_INTEGER,th_block_id,MB_TAG_CREAT|MB_TAG_SPARSE,&def_block
251 const string &name = it->getName();
252 if(name.compare(0,3,
"GEL")!=0)
continue;
254 int block_id = it->getMeshsetId();
256 ostringstream str_order;
257 str_order <<
"block_" << block_id <<
".oRder";
258 config_file_options.add_options()
259 (str_order.str().c_str(),po::value<int>(
260 &block_data[block_id].oRder
261 )->default_value(
order));
263 rval = moab.get_entities_by_handle(
264 it->getMeshset(),block_ents,true
266 material_data[block_id].tEts = block_ents.subset_by_type(MBTET);
267 std::vector<int> block_id_vec(material_data[block_id].tEts.size(),it->getMeshsetId());
269 rval = moab.tag_set_data(
270 th_block_id,material_data[block_id].tEts,&*block_id_vec.begin()
273 ostringstream str_g_alpha;
274 str_g_alpha <<
"block_" << block_id <<
".gAlpha";
275 config_file_options.add_options()
276 (str_g_alpha.str().c_str(),po::value<double>(&material_data[block_id].gAlpha)->default_value(1));
277 ostringstream str_v_alpha;
278 str_v_alpha <<
"block_" << block_id <<
".vAlpha";
279 config_file_options.add_options()
280 (str_v_alpha.str().c_str(),po::value<double>(&material_data[block_id].vAlpha)->default_value(0));
281 ostringstream str_g_beta;
282 str_g_beta <<
"block_" << block_id <<
".gBeta";
283 config_file_options.add_options()
284 (str_g_beta.str().c_str(),po::value<double>(&material_data[block_id].gBeta)->default_value(1));
285 ostringstream str_v_beta;
286 str_v_beta <<
"block_" << block_id <<
".vBeta";
287 config_file_options.add_options()
288 (str_v_beta.str().c_str(),po::value<double>(&material_data[block_id].vBeta)->default_value(0));
289 ostringstream str_g_beta_hat;
290 str_g_beta_hat <<
"block_" << block_id <<
".gBetaHat";
291 config_file_options.add_options()
292 (str_g_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].gBetaHat)->default_value(1));
293 ostringstream str_v_beta_hat;
294 str_v_beta_hat <<
"block_" << block_id <<
".vBetaHat";
295 config_file_options.add_options()
296 (str_v_beta_hat.str().c_str(),po::value<double>(&material_data[block_id].vBetaHat)->default_value(0));
297 ostringstream str_omega;
298 str_omega <<
"block_" << block_id <<
".oMega";
299 config_file_options.add_options()
300 (str_omega.str().c_str(),po::value<double>(&material_data[block_id].oMega)->default_value(1));
301 ostringstream str_viscosity;
302 str_viscosity <<
"block_" << block_id <<
".vIscosity";
303 config_file_options.add_options()
304 (str_viscosity.str().c_str(),po::value<double>(&material_data[block_id].vIscosity)->default_value(1));
305 ostringstream str_permeability;
306 str_permeability <<
"block_" << block_id <<
".pErmeability";
307 config_file_options.add_options()
308 (str_permeability.str().c_str(),po::value<double>(&material_data[block_id].pErmeability)->default_value(1));
309 ostringstream str_mu0;
310 str_mu0 <<
"block_" << block_id <<
".mU0";
311 config_file_options.add_options()
312 (str_mu0.str().c_str(),po::value<double>(&material_data[block_id].mU0)->default_value(1));
316 po::parsed_options parsed = parse_config_file(ini_file,config_file_options,
true);
321 const string &name = it->getName();
322 if(name.compare(0,3,
"GEL")!=0)
continue;
324 if(block_data[it->getMeshsetId()].oRder == -1)
continue;
325 if(block_data[it->getMeshsetId()].oRder ==
order)
continue;
326 PetscPrintf(PETSC_COMM_WORLD,
"Set block %d order to %d\n",it->getMeshsetId(),block_data[it->getMeshsetId()].oRder);
329 Range ents_to_set_order;
330 rval = moab.get_adjacencies(block_ents,3,
false,ents_to_set_order,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
331 ents_to_set_order = ents_to_set_order.subset_by_type(MBTET);
332 rval = moab.get_adjacencies(block_ents,2,
false,ents_to_set_order,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
333 rval = moab.get_adjacencies(block_ents,1,
false,ents_to_set_order,moab::Interface::UNION);
CHKERRQ_MOAB(
rval);
335 ierr = m_field.
set_field_order(ents_to_set_order,
"SPATIAL_POSITION",block_data[it->getMeshsetId()].oRder); CHKERRQ(
ierr);
336 ierr = m_field.
set_field_order(ents_to_set_order,
"CHEMICAL_LOAD",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
337 ierr = m_field.
set_field_order(ents_to_set_order,
"HAT_EPS",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
338 ierr = m_field.
set_field_order(ents_to_set_order,
"SPATIAL_POSITION_DOT",block_data[it->getMeshsetId()].oRder); CHKERRQ(
ierr);
339 ierr = m_field.
set_field_order(ents_to_set_order,
"CHEMICAL_LOAD_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
340 ierr = m_field.
set_field_order(ents_to_set_order,
"HAT_EPS_DOT",block_data[it->getMeshsetId()].oRder-1); CHKERRQ(
ierr);
342 vector<string> additional_parameters;
343 additional_parameters = collect_unrecognized(parsed.options,po::include_positional);
344 for(vector<string>::iterator vit = additional_parameters.begin();
345 vit!=additional_parameters.end();vit++) {
346 ierr = PetscPrintf(PETSC_COMM_WORLD,
"** WARNING Unrecognised option %s\n",vit->c_str()); CHKERRQ(
ierr);
348 }
catch (
const std::exception& ex) {
350 ss << ex.what() << endl;
358 ierr = PetscPrintf(PETSC_COMM_WORLD,
"** Block Id %d\n",mit->first); CHKERRQ(
ierr);
359 PetscPrintf(PETSC_COMM_WORLD,
"vAlpha = %5.4g\n",mit->second.vAlpha);
360 PetscPrintf(PETSC_COMM_WORLD,
"gAlpha = %5.4g\n",mit->second.gAlpha);
361 PetscPrintf(PETSC_COMM_WORLD,
"vBeta = %5.4g\n",mit->second.vBeta);
362 PetscPrintf(PETSC_COMM_WORLD,
"gBeta = %5.4g\n",mit->second.gBeta);
363 PetscPrintf(PETSC_COMM_WORLD,
"vBetaHat = %5.4g\n",mit->second.vBetaHat);
364 PetscPrintf(PETSC_COMM_WORLD,
"gBetaHat = %5.4g\n",mit->second.gBetaHat);
365 PetscPrintf(PETSC_COMM_WORLD,
"pErmeability = %5.4g\n",mit->second.pErmeability);
366 PetscPrintf(PETSC_COMM_WORLD,
"vIscosity = %5.4g\n",mit->second.vIscosity);
367 PetscPrintf(PETSC_COMM_WORLD,
"oMega = %5.4g\n",mit->second.oMega);
368 PetscPrintf(PETSC_COMM_WORLD,
"mU0 = %5.4g\n",mit->second.mU0);
375 if (!check_if_spatial_field_exist) {
377 ierr = m_field.
loop_dofs(
"MESH_NODE_POSITIONS", ent_method_material); CHKERRQ(
ierr);
383 map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
384 mit!=material_data.end();mit++
389 ierr = m_field.set_field(mit->second.mU0,MBVERTEX,vertices,
"CHEMICAL_LOAD"); CHKERRQ(
ierr);
415 map<int,Gel::BlockMaterialData>::iterator mit = material_data.begin();
416 mit!=material_data.end();
423 ierr = MetaNeummanForces::addNeumannBCElements(m_field,
"SPATIAL_POSITION"); CHKERRQ(
ierr);
439 if(it->getName().compare(0,18,
"FLUX_CHEMICAL_LOAD") == 0) {
441 ierr = it->getAttributes(data); CHKERRQ(
ierr);
443 SETERRQ(PETSC_COMM_SELF,1,
"Data inconsistency");
448 strcpy(set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.name,
"HeatFlu");
449 set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.flag1 = 1;
450 set_of_solvent_fluxes[it->getMeshsetId()].dAta.data.value1 = data[0];
453 it->meshset,MBTRI,set_of_solvent_fluxes[it->getMeshsetId()].tRis,
true
456 set_of_solvent_fluxes[it->getMeshsetId()].tRis,MBTRI,
"CHEMICAL_LOAD_FLUX_FE"
472 tags.push_back(Gel::STRESSTOTAL);
473 tags.push_back(Gel::SOLVENTFLUX);
474 tags.push_back(Gel::SOLVENTRATE);
475 tags.push_back(Gel::RESIDUALSTRAINHAT);
488 common_data.spatialPositionNameDot =
"SPATIAL_POSITION_DOT";
489 common_data.muName =
"CHEMICAL_LOAD";
490 common_data.muNameDot =
"CHEMICAL_LOAD_DOT";
491 common_data.strainHatName =
"HAT_EPS";
492 common_data.strainHatNameDot =
"HAT_EPS_DOT";
497 for(
int ss = 0;ss<2;ss++) {
567 DMType dm_name =
"DMGEL";
570 ierr = DMCreate(PETSC_COMM_WORLD,&dm);CHKERRQ(
ierr);
571 ierr = DMSetType(dm,dm_name);CHKERRQ(
ierr);
573 ierr = DMSetFromOptions(dm); CHKERRQ(
ierr);
587 ierr = VecDuplicate(T,&
F); CHKERRQ(
ierr);
593 m_field,
"SPATIAL_POSITION",
A,T,
F
595 spatial_position_bc.methodsOp.push_back(
new TimeForceScale(
"-my_displacements_history",
false));
597 m_field,
"CHEMICAL_LOAD",
"CHEMICAL_LOAD",
A,T,
F
599 concentration_bc.methodsOp.push_back(
new TimeForceScale(
"-my_chemical_load_history",
false));
602 boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
603 boost::ptr_map<string,NodalForce> nodal_forces;
604 boost::ptr_map<string,EdgeForce> edge_forces;
607 ierr = MetaNeummanForces::setMomentumFluxOperators(m_field,neumann_forces,PETSC_NULL,
"SPATIAL_POSITION"); CHKERRQ(
ierr);
613 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
614 mit!=neumann_forces.end();mit++
616 mit->second->methodsOp.push_back(
new TimeForceScale(
"-my_load_history",
false));
619 boost::ptr_map<string,NodalForce>::iterator mit = nodal_forces.begin();
620 mit!=nodal_forces.end();mit++
622 mit->second->methodsOp.push_back(
new TimeForceScale(
"-my_load_history",
false));
625 boost::ptr_map<string,EdgeForce>::iterator mit = edge_forces.begin();
626 mit!=edge_forces.end();mit++
628 mit->second->methodsOp.push_back(
new TimeForceScale(
"-my_load_history",
false));
636 map<int,ThermalElement::FluxData>::iterator sit = set_of_solvent_fluxes.begin();
637 for(;sit!=set_of_solvent_fluxes.end();sit++) {
653 boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
654 for(;mit!=neumann_forces.end();mit++) {
659 boost::ptr_map<string,NodalForce>::iterator fit = nodal_forces.begin();
660 for(;fit!=nodal_forces.end();fit++) {
665 boost::ptr_map<string,EdgeForce>::iterator fit = edge_forces.begin();
666 for(;fit!=edge_forces.end();fit++) {
686 ierr = TSCreate(PETSC_COMM_WORLD,&ts); CHKERRQ(
ierr);
687 ierr = TSSetType(ts,TSBEULER); CHKERRQ(
ierr);
700 ierr = TSSetIFunction(ts,
F,PETSC_NULL,PETSC_NULL); CHKERRQ(
ierr);
701 ierr = TSSetIJacobian(ts,
A,
A,PETSC_NULL,PETSC_NULL); CHKERRQ(
ierr);
710 ts_ctx->get_postProcess_to_do_Monitor().push_back(&post_proc);
711 ierr = TSMonitorSet(ts,f_TSMonitorSet,
ts_ctx,PETSC_NULL); CHKERRQ(
ierr);
715 ierr = TSSetDuration(ts,PETSC_DEFAULT,ftime); CHKERRQ(
ierr);
716 ierr = TSSetFromOptions(ts); CHKERRQ(
ierr);
717 ierr = TSSetDM(ts,dm); CHKERRQ(
ierr);
718 #if PETSC_VERSION_GE(3,7,0)
719 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER); CHKERRQ(
ierr);
721 ierr = TSSolve(ts,T); CHKERRQ(
ierr);
722 ierr = TSGetTime(ts,&ftime); CHKERRQ(
ierr);
723 PetscInt steps,snesfails,rejects,nonlinits,linits;
724 ierr = TSGetTimeStepNumber(ts,&steps); CHKERRQ(
ierr);
725 ierr = TSGetSNESFailures(ts,&snesfails); CHKERRQ(
ierr);
726 ierr = TSGetStepRejections(ts,&rejects); CHKERRQ(
ierr);
727 ierr = TSGetSNESIterations(ts,&nonlinits); CHKERRQ(
ierr);
728 ierr = TSGetKSPIterations(ts,&linits); CHKERRQ(
ierr);
729 PetscPrintf(PETSC_COMM_WORLD,
730 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits %D, linits %D\n",
731 steps,rejects,snesfails,ftime,nonlinits,linits
737 ierr = VecSum(T,&sum); CHKERRQ(
ierr);
738 ierr = PetscPrintf(PETSC_COMM_WORLD,
"sum = %9.8e\n",sum); CHKERRQ(
ierr);
740 ierr = VecNorm(T,NORM_2,&fnorm); CHKERRQ(
ierr);
741 ierr = PetscPrintf(PETSC_COMM_WORLD,
"fnorm = %9.8e\n",fnorm); CHKERRQ(
ierr);
743 if(fabs(sum-7.57417437e+01)>1e-6) {
746 if(fabs(fnorm-3.67577050e+01)>1e-6) {
753 ierr = DMDestroy(&dm); CHKERRQ(
ierr);
754 ierr = VecDestroy(&T); CHKERRQ(
ierr);