Model.hxx 20.8 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
///
////// \file
///
///
/// Created by Sebastien Gilles <sebastien.gilles@inria.fr> on the Wed, 30 Oct 2013 17:29:19 +0100
/// Copyright (c) Inria. All rights reserved.
///
/// \ingroup ModelGroup
/// \addtogroup ModelGroup
/// \{
11

12 13
#ifndef MOREFEM_x_MODEL_x_MODEL_HXX_
# define MOREFEM_x_MODEL_x_MODEL_HXX_
14 15


16
namespace MoReFEM
17
{
18 19


20
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
21
    template<class InputParameterDataT>
22
    Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::Model(const Wrappers::Mpi& a_mpi,
23 24
                                                                 const InputParameterDataT& input_parameter_data,
                                                                 create_domain_list_for_coords a_create_domain_list_for_coords)
25
    : Crtp::CrtpMpi<Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>>(a_mpi)
26
    {
27
        if (a_create_domain_list_for_coords == create_domain_list_for_coords::yes)
28
            Coords::SetCreateDomainListForCoords();
29

30 31 32
        // current date and time system
        time_t now = std::time(0);
        char* date_time = ctime(&now);
33

34
        namespace IPL = Utilities::InputParameterListNS;
35
        using TimeManager = InputParameter::TimeManager;
36
        using Result = InputParameter::Result;
37

38 39
        decltype(auto) time_evolution_policy =
            IPL::Extract<TimeManager::TimeEvolutionPolicy>::Value(input_parameter_data);
40

41 42 43 44 45
        // \todo #574 Add new time evolution policy... To make the code cleaner we would use a factory with policy registering,
        // but as there won't be that much policies it's not that huge a problem to do it manually and edit the if here
        // whenever a new one is added.
        if (time_evolution_policy == "constant_time_step")
            time_manager_ = std::make_unique<TimeManagerInstance<TimeManagerNS::Policy::ConstantTimeStep>>(input_parameter_data);
46 47
        else if (time_evolution_policy == "variable_time_step")
            time_manager_ = std::make_unique<TimeManagerInstance<TimeManagerNS::Policy::VariableTimeStep>>(input_parameter_data);
48 49 50 51
        else
        {
            assert(false && "Unknown time evolution policy!");
        }
52

53 54
        display_value_ = IPL::Extract<Result::DisplayValue>::Value(input_parameter_data);

55
        const auto& mpi = this->GetMpi();
56

57
        output_directory_ = IPL::Extract<Result::OutputDirectory>::Folder(input_parameter_data);
58

59
        if (mpi.IsRootProcessor())
60
        {
61 62
            {
                std::string target(output_directory_);
63
                target += "/input_parameter_file.lua";
64

65 66
                FilesystemNS::File::Copy(input_parameter_data.GetInputFile(),
                                         target,
67
                                         FilesystemNS::File::fail_if_already_exist::no,
68 69
                                         FilesystemNS::File::autocopy::yes,
                                         __FILE__, __LINE__);
70
            }
71

72
            // Specify in output directory how many processors are involved.
73

74 75
            {
                std::string target(output_directory_);
76
                target += "/mpi.hhdata";
77

78
                std::ofstream out;
79
                FilesystemNS::File::Create(out, target, __FILE__, __LINE__);
80

81 82
                out << "Nprocessor: " << mpi.template Nprocessor<int>() << std::endl;
            }
83

84 85 86 87 88
            // Specify the name of the model. It might be useful as same Lua file might be used with different models
            // (e.g. different time scheme or hyperelastic law for hyperelastic model).
            {
                std::string target(output_directory_);
                target += "/model_name.data";
89

90 91
                std::ofstream out;
                FilesystemNS::File::Create(out, target, __FILE__, __LINE__);
92

93 94
                out << DerivedT::ClassName() << std::endl;
            }
95
        }
96

97

98 99
        Wrappers::Petsc::PrintMessageOnFirstProcessor("\n================================================================\n",
                                                      mpi, __FILE__, __LINE__);
100

101
        std::ostringstream oconv;
102
        oconv << "MoReFEM ";
103 104 105 106
        oconv << 8 * sizeof(void*);
        oconv << " bits ";
        Wrappers::Petsc::PrintMessageOnFirstProcessor(oconv.str().c_str(),
                                                      mpi, __FILE__, __LINE__);
107

108 109
        if (mpi.template Nprocessor<int>() == 1)
            Wrappers::Petsc::PrintMessageOnFirstProcessor("on 1 processor \n",
110
                                                          mpi, __FILE__, __LINE__);
111 112 113
        else
            Wrappers::Petsc::PrintMessageOnFirstProcessor("on %d processors \n",
                                                          mpi, __FILE__, __LINE__, mpi.template Nprocessor<int>());
114

115 116
        Wrappers::Petsc::PrintMessageOnFirstProcessor("%s=================================================================\n",
                                                      mpi, __FILE__, __LINE__, date_time);
117
    }
118 119


120 121
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::~Model()
122
    {
123
        const auto& mpi = this->GetMpi();
124

125 126
        Wrappers::Petsc::PrintMessageOnFirstProcessor("\nIf no exception, all results have been printed in %s.\n",
                                                      mpi, __FILE__, __LINE__, GetOutputDirectory().c_str());
127

128 129
        Wrappers::Petsc::PrintMessageOnFirstProcessor("\n==============================================================\n",
                                                      mpi, __FILE__, __LINE__);
130
        Wrappers::Petsc::PrintMessageOnFirstProcessor("MoReFEM %s ended (if exception thrown it will appear afterwards).\n",
131
                                                      mpi, __FILE__, __LINE__, DerivedT::ClassName().c_str());
132

133 134 135
        // current date and time system
        time_t now = time(0);
        char* date_time = ctime(&now);
136

137 138
        Wrappers::Petsc::PrintMessageOnFirstProcessor("%s==============================================================\n",
                                                      mpi, __FILE__, __LINE__, date_time);
139

140
        #ifdef MOREFEM_CHECK_UPDATE_GHOSTS_CALL_RELEVANCE
141
        if (mpi.template Nprocessor<int>() > 1 && mpi.IsRootProcessor())
142
            ::MoReFEM::Internal::Wrappers::Petsc::CheckUpdateGhostManager::GetInstance().Print();
143
        #endif // MOREFEM_CHECK_UPDATE_GHOSTS_CALL_RELEVANCE
144
    }
145 146


147 148
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    bool Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::HasFinished()
149
    {
150
        auto& time_manager = GetNonCstTimeManager();
151

152
        if (time_manager.HasFinished())
153 154
            return true;

155
        if (static_cast<const DerivedT&>(*this).SupplHasFinishedConditions())
156 157 158
            return true;

        return false;
159
    }
160

161

162
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
163
    template<class InputParameterDataT>
164 165 166
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>
    ::Initialize(const InputParameterDataT& input_parameter_data,
                 create_output_dir do_create_output_dir)
167
    {
168
        const auto& mpi = this->GetMpi();
169 170
        decltype(auto) output_directory = GetOutputDirectory();
        const bool is_root_processor = mpi.IsRootProcessor();
171

172
        // Init all the relevant objects found in the input parameter file.
173
        // Ordering is important here: do not modify it lightly!
174

175 176 177 178
        {
            auto& manager = Internal::MeshNS::GeometricMeshRegionManager::CreateOrGetInstance();
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }
179

180
        auto mesh_directory_storage = CreateMeshDataDirectory(do_create_output_dir);
181

182 183
        if (is_root_processor)
            Internal::MeshNS::WriteInterfaceListForEachMesh(mesh_directory_storage);
184

185
        mpi.Barrier();
186

187
        {
188
            auto& manager = DomainManager::CreateOrGetInstance();
189 190
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }
191 192


193
        {
194
            auto& manager = Advanced::LightweightDomainListManager::CreateOrGetInstance();
195 196 197
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }

198
        if (Coords::GetCreateDomainListForCoords() == create_domain_list_for_coords::yes)
199
            CreateDomainListForCoords();
200

201
        // Advanced::SetFromInputParameterData<Internal::PseudoNormalsManager>(input_parameter_data); #938 not activated for the moment.
202 203 204 205 206 207 208 209 210 211
        {
            auto& manager = UnknownManager::CreateOrGetInstance();
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }

        {
            auto& manager = DirichletBoundaryConditionManager::CreateOrGetInstance();
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }

212

213 214 215 216
        // Write the list of unknowns.
        // #289 We assume here there is one model, as the target file gets a hardcoded name...
        if (is_root_processor)
            WriteUnknownList(output_directory);
217

218 219 220 221 222 223 224 225 226 227
        {
            auto& manager = Internal::NumberingSubsetNS::NumberingSubsetManager::CreateOrGetInstance();
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }

        {
            auto& manager = GodOfDofManager::CreateOrGetInstance();
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager, mpi);
        }

228

229 230
        Internal::ModelNS::InitGodOfDof(input_parameter_data,
                                        DoConsiderProcessorWiseLocal2GlobalT,
231
                                        std::move(mesh_directory_storage));
232

233 234
        // As FiberListManager gets a constructor with arguments, call it explicitly there.
        decltype(auto) time_manager = GetTimeManager();
235

236 237 238 239 240 241 242 243 244 245
        {
            auto& manager = Internal::FiberNS::FiberListManager<ParameterNS::Type::scalar>::CreateOrGetInstance(time_manager);
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }

        {
            auto& manager = Internal::FiberNS::FiberListManager<ParameterNS::Type::vector>::CreateOrGetInstance(time_manager);
            Advanced::SetFromInputParameterData<>(input_parameter_data, manager);
        }

246

247
        static_cast<DerivedT&>(*this).SupplInitialize(input_parameter_data);
248

249 250
        if (do_clear_god_of_dof_temporary_data_after_initialize_)
            ClearGodOfDofTemporaryData();
251

252
        ClearAllBoundaryConditionInitialValueList();
253
    }
254 255


256 257
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::InitializeStep()
258
    {
259 260
        // Update time for current time step.
        UpdateTime();
261

262
        // Print time information
263
        if (do_print_new_time_iteration_banner_)
264
        {
265 266 267 268
            if ((this->GetTimeManager().NtimeModified() % GetDisplayValue()) == 0)
            {
                PrintNewTimeIterationBanner();
            }
269
        }
270

271
        // Additional operations required by DerivedT.
272
        static_cast<DerivedT&>(*this).SupplInitializeStep();
273
    }
274 275


276 277
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::FinalizeStep()
278 279
    {
        // Additional operations required by DerivedT.
280
        static_cast<DerivedT&>(*this).SupplFinalizeStep();
281
    }
282 283


284 285 286
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::Finalize()
    {
287
        const auto& mpi = this->GetMpi();
288

289
        decltype(auto) time_keep = TimeKeep::GetInstance();
290

291
        std::string time_end = time_keep.TimeElapsedSinceBeginning();
292

293 294 295 296 297 298 299
        Wrappers::Petsc::PrintMessageOnFirstProcessor("\n----------------------------------------------\n",
                                                      mpi, __FILE__, __LINE__);
        Wrappers::Petsc::PrintMessageOnFirstProcessor("Time of execution : %s.\n",
                                                      mpi, __FILE__, __LINE__,
                                                      time_end.c_str());
        Wrappers::Petsc::PrintMessageOnFirstProcessor("----------------------------------------------\n",
                                                      mpi, __FILE__, __LINE__);
300

301 302 303
        // Destroy manually some singletons that would blow up at the end of the program: due to a dependancy
        // their content must be destroyed before the call to Mpi::Finalize(), which occurs before the natural end of
        // singletons.
304 305
        Internal::FiberNS::FiberListManager<ParameterNS::Type::scalar>::GetInstance().Destroy();
        Internal::FiberNS::FiberListManager<ParameterNS::Type::vector>::GetInstance().Destroy();
306

307
        time_keep.PrintEndProgram();
308

309 310 311 312
        // Additional operations required by DerivedT.
        static_cast<DerivedT&>(*this).SupplFinalize();
    }

313

314
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
315
    template<class InputParameterDataT>
316 317
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::Run(const InputParameterDataT& input_parameter_data,
                                                                    create_output_dir do_create_output_dir)
318
    {
319
        auto& crtp_helper = static_cast<DerivedT&>(*this);
320

321
        crtp_helper.Initialize(input_parameter_data, do_create_output_dir);
322

323
        while (!crtp_helper.HasFinished())
324 325 326 327 328
        {
            crtp_helper.InitializeStep();
            crtp_helper.Forward();
            crtp_helper.FinalizeStep();
        }
329

330 331
        crtp_helper.Finalize();
    }
332 333


334 335 336
    ////////////////////////
    // ACCESSORS          //
    ////////////////////////
337 338 339



340 341
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    inline const GeometricMeshRegion& Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>
342
    ::GetGeometricMeshRegion(unsigned int mesh_index) const
343
    {
344
        return Internal::MeshNS::GeometricMeshRegionManager::GetInstance().GetMesh(mesh_index);
345
    }
346 347


348 349
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    inline GeometricMeshRegion& Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>
350
    ::GetNonCstGeometricMeshRegion(unsigned int mesh_index) const
351
    {
352
        return Internal::MeshNS::GeometricMeshRegionManager::GetInstance().GetNonCstMesh(mesh_index);
353
    }
354

355

356
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
GILLES Sebastien's avatar
GILLES Sebastien committed
357
    const GodOfDof& Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::GetGodOfDof(unsigned int unique_id) const
358
    {
GILLES Sebastien's avatar
GILLES Sebastien committed
359
        return GodOfDofManager::GetInstance().GetGodOfDof(unique_id);
360
    }
361

362

363 364
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    inline TimeManager& Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::GetNonCstTimeManager()
365
    {
366 367
        assert(!(!time_manager_));
        return *time_manager_;
368
    }
369 370


371 372
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    inline const TimeManager& Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::GetTimeManager() const
373
    {
374 375
        assert(!(!time_manager_));
        return *time_manager_;
376
    }
377 378 379



380 381
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    const std::string& Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::GetOutputDirectory() const
382 383 384 385 386
    {
        assert(!output_directory_.empty());
        return output_directory_;
    }

387 388


389 390 391
    ////////////////////////
    // PRIVATE METHODS    //
    ////////////////////////
392 393


394

395 396
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::UpdateTime()
397
    {
398
        time_manager_->IncrementTime();
399
    }
400 401


402 403
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::PrintNewTimeIterationBanner() const
404
    {
405
        const auto& mpi = this->GetMpi();
406

407
        Wrappers::Petsc::PrintMessageOnFirstProcessor("\n---------------------------------------------------\n",
408
                                                      mpi, __FILE__, __LINE__);
409
        Wrappers::Petsc::PrintMessageOnFirstProcessor("%s Iteration: %d, Time: %f -> %f \n",
410
                                                      mpi, __FILE__, __LINE__,
411
                                                      DerivedT::ClassName().c_str(),
412
                                                      time_manager_->NtimeModified(),
413
                                                      time_manager_->GetTime() - time_manager_->GetTimeStep(),
414
                                                      time_manager_->GetTime());
415
        Wrappers::Petsc::PrintMessageOnFirstProcessor("---------------------------------------------------\n",
416 417
                                                      mpi, __FILE__, __LINE__);
    }
418 419


420 421 422 423 424
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    unsigned int Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::GetDisplayValue() const
    {
        return display_value_;
    }
425 426


427 428 429 430 431 432 433 434
    template
    <
        class DerivedT,
        DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT
    >
    std::map<unsigned int, std::string> Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>
    ::CreateMeshDataDirectory(create_output_dir do_create_output_dir) const
    {
435

436
        const auto& mpi = this->GetMpi();
437

438
        std::map<unsigned int, std::string> ret;
439

440 441 442 443
        decltype(auto) mesh_manager = Internal::MeshNS::GeometricMeshRegionManager::GetInstance();

        decltype(auto) mesh_id_list = mesh_manager.GetUniqueIdList();
        decltype(auto) output_dir = GetOutputDirectory();
444

445 446
        for (const auto mesh_id : mesh_id_list)
            ret.insert({ mesh_id, output_dir + "/Mesh_" + std::to_string(mesh_id)});
447

448 449 450 451 452
        if (do_create_output_dir == create_output_dir::yes && mpi.IsRootProcessor())
        {
            for (const auto& pair : ret)
            {
                const auto& mesh_directory = pair.second;
453

454 455
                if (FilesystemNS::Folder::DoExist(mesh_directory))
                    FilesystemNS::Folder::Remove(mesh_directory, __FILE__, __LINE__);
456

457 458 459
                FilesystemNS::Folder::Create(mesh_directory, __FILE__, __LINE__);
            }
        }
460

461
        mpi.Barrier();
462

463 464
        return ret;
    }
465 466


467 468 469 470 471
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::SetClearGodOfDofTemporaryDataToFalse()
    {
        do_clear_god_of_dof_temporary_data_after_initialize_ = false;
    }
GILLES Sebastien's avatar
GILLES Sebastien committed
472

473

474 475 476 477 478
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    inline void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::SetDoPrintNewTimeIterationBanner(bool do_print) noexcept
    {
        do_print_new_time_iteration_banner_ = do_print;
    }
479 480


481 482 483 484 485
    template<class DerivedT, DoConsiderProcessorWiseLocal2Global DoConsiderProcessorWiseLocal2GlobalT>
    inline void Model<DerivedT, DoConsiderProcessorWiseLocal2GlobalT>::CreateDomainListForCoords()
    {
        {
            auto& domain_manager = DomainManager::CreateOrGetInstance();
486

487
            const auto& domain_list = domain_manager.GetDomainList();
488

489 490 491
            for (const auto& domain_and_id : domain_list)
            {
                const auto& domain_ptr = domain_and_id.second;
492

493
                assert(!(!domain_ptr));
494

495
                const auto& domain = *domain_ptr;
496

497
                const auto domain_unique_id = domain_and_id.first;
498

499
                const auto& mesh = domain.GetGeometricMeshRegion();
500

501
                const auto& geometric_elt_list = mesh.GetGeometricEltList();
502

503 504 505
                for (const auto& geometric_elt_ptr : geometric_elt_list)
                {
                    assert(!(!geometric_elt_ptr));
506

507
                    const auto& geom_elt = *geometric_elt_ptr;
508

509 510 511
                    if (domain.IsGeometricEltInside(geom_elt))
                    {
                        auto& coords_list_elem = geom_elt.GetCoordsList();
512

513 514 515
                        for (auto& coords_ptr : coords_list_elem)
                        {
                            assert(!(!coords_ptr));
516

517
                            auto& coords = *coords_ptr;
518

519 520 521 522 523 524 525
                            coords.AddDomainContainingCoords(domain_unique_id);
                        }
                    }
                }
            }
        }
    }
526 527


528
} // namespace MoReFEM
529

530

531 532 533
/// @} // addtogroup ModelGroup


534
#endif // MOREFEM_x_MODEL_x_MODEL_HXX_