Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Subdomain deactivation within time intervals #2297

Merged
merged 27 commits into from Jan 8, 2019

Conversation

Projects
None yet
3 participants
@wenqing
Copy link
Member

wenqing commented Dec 13, 2018

This PR presents an important feature of subdomain deactivation within time intervals. As a benefit of the feature, the modelling of time dependent excavation is possible.

Based on PR #2244, this PR extends the subdomain deactivation as a general approach for all processes and removes the linear solver restriction. All comments and suggestions to PR #2244 were considered. The implementation includes:

  1. subdomain deactivation data are set as members of ProcessVariable.
  2. nodes in the deactivated subdomains, which are excluded those on the borders among deactivated and activated subdomains, are assigned with zero Dirichlet BC in order to guarantee that the system of equations is not singular. For this, a constant zero Parameter is created.
  3. a boolean vector (belong to ProcessVariable) is used to record the element status of deactivation. The vector is only initialized once for each time interval.
  4. two new members, executeSelectedMemberDereferenced and executeMemberOnDereferenced,
    were added to struct SerialExecutor to use the boolean vector of ProcessVariable to skip the elements that are deactivated.
  5. In all process classes, executeMemberDereferenced and executeOnDereferenced are replaced with executeSelectedMemberDereferenced and executeMemberOnDereferenced, respectively.

So far, only SmallDeformation gets a benchmark for this feature in this PR. The benchmark is the same as that in PR #2244, a comparison of the results of the present benchmark (left) and disc_with_hole (right):
47353534-ca3d2280-d6bc-11e8-90fb-5838b9e42af8

Benchmarks for other processes with this feature will be presented later on.

@wenqing wenqing force-pushed the wenqing:new_deactivate_doms branch from df3fac5 to c3fa725 Dec 13, 2018

@@ -30,6 +30,12 @@ class TimeInterval final
{
}

explicit TimeInterval(const TimeInterval& time_inverval)

This comment has been minimized.

@endJunction

endJunction Dec 13, 2018

Member

explicit is definitely not needed, maybe = default; would be sufficient too. operator=() missing.

This comment has been minimized.

@wenqing

wenqing Dec 13, 2018

Author Member

took = default. Why need the assignment operator?

This comment has been minimized.

@wenqing

wenqing Dec 18, 2018

Author Member

added a default assignment operator.

Method method, Container const& container,
std::vector<bool> const& container_selector, Args&&... args)
{
if(container_selector.empty())

This comment has been minimized.

@endJunction

endJunction Dec 13, 2018

Member

For discussion: I'd rather assert that the container_selector is not empty.

@@ -34,20 +34,44 @@ DirichletBoundaryConditionWithinTimeInterval::
int const variable_id, int const component_id)
: _parameter(parameter),
_bc_mesh(bc_mesh),
_nodes_in_bc_mesh(bc_mesh.getNodes()),

This comment has been minimized.

@endJunction

endJunction Dec 13, 2018

Member

For discussion: Why is this needed.

This comment has been minimized.

@wenqing

wenqing Dec 18, 2018

Author Member

As explained.

//! \ogs_file_param{prj__process_variables__process_variable__deactivated_subdomains__deactivated_subdomain__material_ids}
config.getConfigParameter<std::vector<int>>("material_ids",
std::vector<int>{});
if (!deactivated_subdomain_material_ids.empty())

This comment has been minimized.

@endJunction

endJunction Dec 13, 2018

Member

Invert the if condition => reduced indentation

This comment has been minimized.

@wenqing

wenqing Dec 13, 2018

Author Member

done.

@@ -24,6 +24,8 @@

#include "BaseLib/Algorithm.h"
#include "BaseLib/FileTools.h"
// include due to the forward declaration in ProcessVariable.h
#include "BaseLib/TimeInterval.h"

This comment has been minimized.

@endJunction

endJunction Dec 13, 2018

Member

This looks to me to be at wrong place. It should be in DeactivatedSubdomain.h because of the unique_ptr used there; see https://stackoverflow.com/a/6089065

This comment has been minimized.

@wenqing

wenqing Dec 13, 2018

Author Member

took it.

@wenqing wenqing force-pushed the wenqing:new_deactivate_doms branch from c3fa725 to f70d099 Dec 13, 2018

std::vector<std::unique_ptr<DeactivatedSubdomain const>>
_deactivated_subdomains;

mutable std::vector<bool> _element_deactivation_flags;

This comment has been minimized.

@wenqing

wenqing Dec 14, 2018

Author Member

Will change this to a vector of activated element IDs to remove if-condition inside the loop of the executor.

This comment has been minimized.

@wenqing

wenqing Dec 14, 2018

Author Member

Since it is a simple change, I will make it after you finish reviewing the present commits.

This comment has been minimized.

@wenqing

wenqing Dec 18, 2018

Author Member

Changed for a test (see commit f6afbe4). A little acceleration in assembly loop but memory increases 64 times (on 64 bit machine) with the replacement of a boolean vector with a std::size_t vector. See the efficiency comparison below.

@wenqing

This comment has been minimized.

Copy link
Member Author

wenqing commented Dec 14, 2018

@bilke What is the following error in git-check?

[MSVC] -> found 428 issues (skipped 46 duplicates)
[MSVC] [ERROR] Can't resolve absolute paths for 417 files:
[MSVC] [ERROR] - ../ThirdParty/metis/GKlib/util.c
...

@wenqing wenqing force-pushed the wenqing:new_deactivate_doms branch from f70d099 to 80adcec Dec 17, 2018

@@ -37,16 +36,16 @@ struct SerialExecutor
/// \tparam C input container type.
/// \tparam Args types additional arguments passed to \c f.
///
/// \param f a function that accepts a a reference to container's elements,
/// \param f a function that accepts a a reference to container's

This comment has been minimized.

@endJunction

endJunction Dec 18, 2018

Member

double "a" before reference.

This comment has been minimized.

@wenqing

wenqing Dec 18, 2018

Author Member

removed by amending the last commit.

@wenqing

This comment has been minimized.

Copy link
Member Author

wenqing commented Dec 18, 2018

Comparison of different members of struct SerialExecutor

Test cases

  1. original one, SerialExecutor::executeMemberDereferenced.

  2. new member, SerialExecutor::executeSelectedMemberOnDereferenced, which uses std::vector<bool> container_selector to distinguish which container members have to been skipped. The member contains

       assert(container.size() == container_selector.size());
       for (std::size_t i = 0; i < container.size(); i++)
       {
          if (container_selector[i])
              continue;
          (object.*method)(i, *container[i], std::forward<Args>(args)...);
       }
    

    For the test, all entries of container_selector are set to false.

  3. The same new member as 2. However, the argument of std::vector container_selector was replaced with std::vector<std::size_t> const& active_container_ids, which takes ProcessVariable::_ids_of_active_elements in computations (the last commit: 8610977). The content of executeSelectedMemberOnDereferenced was changed to

       for (std::size_t i = 0; i < active_container_ids.size(); i++)
       {
           (object.*method)(i, *container[active_container_ids[i]],
                            std::forward<Args>(args)...);
       }
    

    which drops the if-condition inside the loop that exists in the implementation in 2.
    For the tests, active_container_ids is set to contain all element IDs.

Benchmark
The assembly part of Mechanics/Linear/square_1e5.prj.

Run time comparison
The comparison uses the release mode compiled exactable files.
| Case 1 | Case 2 | Case 3 |
|---|---|---|
| 0.383908 s | 0.382762 s | 0.382998 s |

One can see that there is no much differences in the runtime of the three cases.

Additional comparison
The tests in the above comparison do not have deactivated subdomains. In this comparison, the benchmark for this PR, which contains deactivated elements, is used. Therefore only case 2 and case 3 are compared. The run times of the entire computations by Case 2 and Case 3 are:
| Case 2 | Case 3 |
|---|---|
| 0.48 s | 0.37 s |
Relatively the computation time with Case 3 decreases (0.48-0.37)/0.48=0.22916666666666666 due to the condensed assembly loop in Case 3. The cost is that the memory for the element deactivation increases 2768 * 64 times of that by Case 2 (2768 * 1 bit of boolean) , where 2768 is the number of the elements in the benchmark.

Conclusions
1 If there is not any deactivated subdomain in computation, the changes of SerialExecutor in this PR does not influence the computational efficiency distinctively.
2. If there are deactivated subdomains in computation, Case 3 shows a significant run time reduction compared to that by Case 2 but with a little bit memory usage increase.

@wenqing wenqing force-pushed the wenqing:new_deactivate_doms branch 4 times, most recently from 6cfaad6 to 8610977 Dec 18, 2018


for (std::size_t i = 0; i < mesh.getNumberOfElements(); i++)
{
if (ids != (*material_ids)[i])

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

Please add an OGS_FATAL if the material_ids are needed but not available. The pointer might be null.

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

Moved the checking of material_ids from the ctor of Process to this function.

continue;

auto* element = mesh.getElement(i);
deactivated_elements.push_back(

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

Would you please have a look if deactivated_elements could be of std::vector<MeshLib::Element const*> type?
Same question to the deactivated_nodes.

Maybe we need a const overload for MeshLib::cloneElements.

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

Tried but have to keep as they are. Reasons:

  1. std::vector<MeshLib::Elemen *>: MeshLib::createMeshFromElementSelection needs this type
  2. std::vector<MeshLib::Node*>: bc_mesh.getNodes() returns this type of variable, which is needed to initialize DeactivatedSubdomain:: std::vector<MeshLib::Node*> const inactive_nodes

If const is added, several more files have to to changed as well. This can be done in another PR.

@@ -568,13 +574,17 @@ void HydroMechanicsProcess<GlobalDim>::assembleWithJacobianConcreteProcess(
{
DBUG("AssembleWithJacobian HydroMechanicsProcess.");

const int process_id =
_use_monolithic_scheme ? 0 : _coupled_solutions->process_id;

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

For another PR; I'd suggest to extract a function currentProcessId() into the base Process class.

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

TODO once this PR is merged.

@@ -568,13 +574,17 @@ void HydroMechanicsProcess<GlobalDim>::assembleWithJacobianConcreteProcess(
{
DBUG("AssembleWithJacobian HydroMechanicsProcess.");

const int process_id =
_use_monolithic_scheme ? 0 : _coupled_solutions->process_id;
ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

For another PR; Maybe also a processVariable(int) function too with appropriate safety checks.

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

TODO once this PR is merged too.

auto const& deactivated_sudomain_meshes =
(*deactivated_subdomain).deactivated_sudomain_meshes;
for (auto const& deactivated_sudomain_mesh :
deactivated_sudomain_meshes)

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

typo; subdomain

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

corrected.

// Copy the time interval.
std::unique_ptr<BaseLib::TimeInterval> time_interval =
std::make_unique<BaseLib::TimeInterval>(
(*(*deactivated_subdomain).time_interval));

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

I'd think that (*(*deactivated_subdomain).time_interval) is equivalent to *deactivated_subdomain->time_interval.
Same pattern below when using (*deactivated_sudomain_mesh).

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

Changed.

This comment has been minimized.

@wenqing

wenqing Jan 8, 2019

Author Member

Changed the remaining stuffs by amending commit b4200a7.

}
}

void ProcessVariable::checkElementDeactivation(double const time)

This comment has been minimized.

@endJunction

endJunction Dec 19, 2018

Member

It was confusing me above in the for-loop where this check is called; The function does not only check but update some internal variables. Then it's better to call it updateDeactivatedSubdomains or similar.

This comment has been minimized.

@wenqing

wenqing Jan 7, 2019

Author Member

Changed.

@endJunction
Copy link
Member

endJunction left a comment

Minor improvements needed.
Additional question: is it working with petsc? If not, add a TODO in the code, please.

@@ -9,14 +9,23 @@

#include "ProcessVariable.h"

#include <iostream>

This comment has been minimized.

@TomFischer

TomFischer Jan 7, 2019

Member

Is iostream needed?

This comment has been minimized.

@wenqing

wenqing Jan 8, 2019

Author Member

removed.

@TomFischer

This comment has been minimized.

Copy link
Member

TomFischer commented Jan 7, 2019

You substituted almost all occurrences of executeMemberOnDereferenced with executeSelectedMemberOnDereferenced. I searched for executeMemberOnDereferenced and found it still in:

  • ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
  • ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
  • ProcessLib/PhaseField/PhaseFieldProcess.cpp
  • ProcessLib/BoundaryCondition/NormalTractionBoundaryCondition-impl.h
  • ProcessLib/BoundaryCondition/GenericNaturalBoundaryCondition-impl.h
  • ProcessLib/BoundaryCondition/Python/PythonBoundaryCondition.cpp
  • ProcessLib/BoundaryCondition/GenericNonuniformNaturalBoundaryCondition-impl.h:

Why is the substitution not necessary in these cases?

@TomFischer
Copy link
Member

TomFischer left a comment

If questions are answered and iostream header removed:

@wenqing wenqing force-pushed the wenqing:new_deactivate_doms branch from 8610977 to a997a6f Jan 7, 2019

@wenqing

This comment has been minimized.

Copy link
Member Author

wenqing commented Jan 7, 2019

@endJunction This feature does not work under PETSc so far. Added TODO in the source code.

@TomFischer Made the replacement in the remaining member of PhaseField, and all associated members of SmallDeformationNonlocal. For the mentioned files about source terms and boundary conditions, the replacement is skipped because the nodes in the deactivated nodes are assigned with zero Dirichlet BC.

wenqing added some commits Dec 3, 2018

[include] added a constant zero-value parameter for the element deact…
…ivation approach.

also added two includes to ogs.cpp and ProjectData.cpp due to the forward declaration of the related classes.
[deactivation] removed the forward declaration of TimeInterval and De…
…activatedSubdomain

due to its std::unique_ptr type instance.
[deactivation] In ProcessVariable, replaced std::vector<bool> _elemen…
…t_deactivation_flags with

std::vector<std::size_t> _ids_of_active_elements.

@wenqing wenqing force-pushed the wenqing:new_deactivate_doms branch from a997a6f to 53bbc5b Jan 8, 2019

@endJunction endJunction merged commit 8bd8e62 into ufz:master Jan 8, 2019

2 of 3 checks passed

continuous-integration/jenkins/pr-merge This commit cannot be built
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
deploy/netlify Deploy preview ready!
Details
/// @param balance The vector the integration results will be stored in.
/// @param t The balance will be computed at the time t.
/// @param bulk_mesh Stores a reference to the bulk mesh that is needed to
/// fetch the information for the integration over the surface mesh.
/// @active_element_ids The IDs of active elements. It is empty if there is

This comment has been minimized.

@endJunction

endJunction Jan 8, 2019

Member

@wenqing Please fix the doxygen comment.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.