[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Fri, 3 Apr 2020 20:09:02 -0400 (EDT) |
branch: master
commit 6ddb89892316a0c6a0cb7eda9a6ab0db07c48ab4
Author: Konstantinos Poulios <address@hidden>
AuthorDate: Wed Jan 29 16:49:18 2020 +0100
Add support for internal variable condensation to the model object
---
src/getfem/getfem_accumulated_distro.h | 6 +-
src/getfem/getfem_models.h | 42 ++++++--
src/getfem_generic_assembly_workspace.cc | 4 +-
src/getfem_model_solvers.cc | 90 +++++++++++++++-
src/getfem_models.cc | 178 +++++++++++++++++++++++--------
5 files changed, 261 insertions(+), 59 deletions(-)
diff --git a/src/getfem/getfem_accumulated_distro.h
b/src/getfem/getfem_accumulated_distro.h
index ebf8226..a428202 100644
--- a/src/getfem/getfem_accumulated_distro.h
+++ b/src/getfem/getfem_accumulated_distro.h
@@ -173,12 +173,16 @@ namespace detail {
}
}
- operator T&(){
+ T& get(){
if (distributed.num_threads() == 1 ||
distributed.this_thread() == 0) return original;
else return distributed;
}
+ operator T&(){ // implicit conversion
+ return get();
+ }
+
T& operator = (const T &x){
return distributed = x;
}
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index c28749f..07f056c 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -121,9 +121,14 @@ namespace getfem {
bool is_linear_;
bool is_symmetric_;
bool is_coercive_;
- mutable model_real_sparse_matrix rTM; // tangent matrix, real version
+ mutable model_real_sparse_matrix
+ rTM, // tangent matrix (only primary variables), real version
+ internal_rTM; // coupling matrix between internal and primary vars (no
empty rows)
mutable model_complex_sparse_matrix cTM; // tangent matrix, complex version
- mutable model_real_plain_vector rrhs;
+ mutable model_real_plain_vector
+ rrhs, // residual vector of primary variables (after
condensation)
+ full_rrhs, // residual vector of primary and internal variables
(pre-condensation)
+ internal_sol; // partial solution for internal variables (after
condensation)
mutable model_complex_plain_vector crhs;
mutable bool act_size_to_be_done;
dim_type leading_dim;
@@ -278,6 +283,10 @@ namespace getfem {
BUILD_ON_DATA_CHANGE = 4,
BUILD_WITH_LIN = 8, // forced calculation of linear terms
BUILD_RHS_WITH_LIN = 9, // = BUILD_RHS | BUILD_WITH_LIN_RHS
+ BUILD_WITH_INTERNAL = 16,
+ BUILD_RHS_WITH_INTERNAL = 17, // = BUILD_RHS | BUILD_WITH_INTERNAL
+ BUILD_MATRIX_CONDENSED = 18, // = BUILD_MATRIX | BUILD_WITH_INTERNAL
+ BUILD_ALL_CONDENSED = 19, // = BUILD_ALL | BUILD_WITH_INTERNAL
};
protected:
@@ -571,7 +580,13 @@ namespace getfem {
bool is_linear() const { return is_linear_; }
/** Total number of degrees of freedom in the model. */
- size_type nb_dof() const;
+ size_type nb_dof(bool with_internal=false) const;
+
+ /** Number of internal degrees of freedom in the model. */
+ size_type nb_internal_dof() const { return nb_dof(true)-nb_dof(); }
+
+ /** Number of primary degrees of freedom in the model. */
+ size_type nb_primary_dof() const { return nb_dof(); }
/** Leading dimension of the meshes used in the model. */
dim_type leading_dimension() const { return leading_dim; }
@@ -898,10 +913,11 @@ namespace getfem {
size_type qdim_of_variable(const std::string &name) const;
/** Gives the access to the tangent matrix. For the real version. */
- const model_real_sparse_matrix &real_tangent_matrix() const {
+ const model_real_sparse_matrix &
+ real_tangent_matrix(bool internal=false) const {
GMM_ASSERT1(!complex_version, "This model is a complex one");
context_check(); if (act_size_to_be_done) actualize_sizes();
- return rTM;
+ return internal ? internal_rTM : rTM;
}
/** Gives the access to the tangent matrix. For the complex version. */
@@ -913,19 +929,27 @@ namespace getfem {
/** Gives access to the right hand side of the tangent linear system.
For the real version. An assembly of the rhs has to be done first. */
- const model_real_plain_vector &real_rhs() const {
+ const model_real_plain_vector &real_rhs(bool with_internal=false) const {
GMM_ASSERT1(!complex_version, "This model is a complex one");
context_check(); if (act_size_to_be_done) actualize_sizes();
- return rrhs;
+ return (with_internal && gmm::vect_size(full_rrhs)) ? full_rrhs : rrhs;
}
/** Gives write access to the right hand side of the tangent linear system.
Some solvers need to manipulate the model rhs directly so that for
example internal condensed variables can be treated properly. */
- model_real_plain_vector &set_real_rhs() const {
+ model_real_plain_vector &set_real_rhs(bool with_internal=false) const {
+ GMM_ASSERT1(!complex_version, "This model is a complex one");
+ context_check(); if (act_size_to_be_done) actualize_sizes();
+ return (with_internal && gmm::vect_size(full_rrhs)) ? full_rrhs : rrhs;
+ }
+
+ /** Gives access to the partial solution for condensed internal variables.
+ A matrix assembly with condensation has to be done first. */
+ const model_real_plain_vector &internal_solution() const {
GMM_ASSERT1(!complex_version, "This model is a complex one");
context_check(); if (act_size_to_be_done) actualize_sizes();
- return rrhs;
+ return internal_sol;
}
/** Gives access to the part of the right hand side of a term of
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index 64183e5..59f2314 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -1001,8 +1001,8 @@ namespace getfem {
nb_tmp_dof(0), macro_dict(md_.macro_dictionary())
{
init();
- nb_prim_dof = with_parent_variables ? md->nb_dof() : 0;
- nb_intern_dof = 0;
+ nb_prim_dof = with_parent_variables ? md->nb_primary_dof() : 0;
+ nb_intern_dof = with_parent_variables ? md->nb_internal_dof() : 0;
if (var_inherit == inherit::ALL) { // enable model's disabled variables
model::varnamelist vlmd;
md->variable_list(vlmd);
diff --git a/src/getfem_model_solvers.cc b/src/getfem_model_solvers.cc
index 0a5c10d..bcaf1af 100644
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@ -191,6 +191,89 @@ namespace getfem {
+ /* ***************************************************************** */
+ /* Non-linear model problem with internal variable condensation. */
+ /* ***************************************************************** */
+ template <typename PLSOLVER>
+ class nonlin_condensed_model_pb : public nonlin_model_pb<PLSOLVER> {
+ public:
+ using typename pb_base<PLSOLVER>::MATRIX;
+ using typename pb_base<PLSOLVER>::VECTOR;
+ using typename pb_base<PLSOLVER>::R;
+
+ private:
+ using nonlin_model_pb<PLSOLVER>::md;
+ bool condensed_vars;
+ const MATRIX &internal_K;
+ VECTOR &full_rhs;
+
+ public:
+ virtual const VECTOR &residual() const { return full_rhs; }
+ // A norm1 seems to be better than a norm2 (at least for contact problems).
+ virtual R residual_norm() { return gmm::vect_norm1(full_rhs); }
+
+ using pb_base<PLSOLVER>::state_vector;
+ using pb_base<PLSOLVER>::state_norm;
+
+ virtual void add_to_residual(VECTOR &extra_rhs, R mult=1.) {
+ if (mult == R(1)) gmm::add(extra_rhs, full_rhs);
+ else if (mult != R(0)) gmm::add(gmm::scaled(extra_rhs, mult), full_rhs);
+ }
+
+ using pb_base<PLSOLVER>::perturbation;
+
+ virtual void linear_solve(VECTOR &dr, gmm::iteration &iter) {
+ VECTOR dr0(md.nb_primary_dof(), R(0));
+ pb_base<PLSOLVER>::linear_solve(dr0, iter);
+ gmm::sub_interval I_prim(0, md.nb_primary_dof()),
+ I_intern(md.nb_primary_dof(), md.nb_internal_dof());
+ gmm::copy(dr0, gmm::sub_vector(dr, I_prim));
+ // recover solution of condensed variables: intern_K*prim_sol +
intern_sol
+ gmm::mult(internal_K,
+ gmm::scaled(gmm::sub_vector(dr, I_prim), scalar_type(-1)),
+ md.internal_solution(),
+ gmm::sub_vector(dr, I_intern));
+ }
+
+ virtual void compute_tangent_matrix() {
+ md.to_variables(state_vector(), condensed_vars);
+ md.assembly(condensed_vars ? model::BUILD_MATRIX_CONDENSED
+ : model::BUILD_MATRIX);
+ }
+
+ virtual void compute_residual() {
+ md.to_variables(state_vector(), condensed_vars);
+ md.assembly(condensed_vars ? model::BUILD_RHS_WITH_INTERNAL // -->
full_rhs
+ : model::BUILD_RHS); // ---> rhs
(= full_rhs)
+ }
+
+ nonlin_condensed_model_pb(model &md_, abstract_newton_line_search &ls_,
+ PLSOLVER linear_solver_) = delete;
+ };
+
+ template <>
+ nonlin_condensed_model_pb<rmodel_plsolver_type>::nonlin_condensed_model_pb
+ (model &md_, abstract_newton_line_search &ls_, rmodel_plsolver_type
linsolv)
+ : nonlin_model_pb<rmodel_plsolver_type>(md_, ls_, linsolv),
+ condensed_vars(md_.has_internal_variables()),
+ internal_K(md_.real_tangent_matrix(condensed_vars)),
+ full_rhs(md_.set_real_rhs(condensed_vars))
+ {
+ gmm::resize(state_vector(), md.nb_dof(condensed_vars));
+ md.from_variables(state_vector(), condensed_vars);
+ }
+
+ template <>
+ nonlin_condensed_model_pb<cmodel_plsolver_type>::nonlin_condensed_model_pb
+ (model &md_, abstract_newton_line_search &ls_, cmodel_plsolver_type
linsolv)
+ : nonlin_model_pb<cmodel_plsolver_type>(md_, ls_, linsolv),
+ condensed_vars(false), internal_K(md_.complex_tangent_matrix()),
+ full_rhs(md_.set_complex_rhs())
+ {
+ GMM_ASSERT1(false, "No support for internal variable condensation in "
+ "complex valued models.");
+ }
+
/* ***************************************************************** */
@@ -311,8 +394,11 @@ namespace getfem {
mdpb.linear_solve(mdpb.state_vector(), iter);
md.to_variables(mdpb.state_vector()); // copy the state vector into the
model variables
} else {
- std::unique_ptr<nonlin_model_pb<PLSOLVER>> mdpb
- = std::make_unique<nonlin_model_pb<PLSOLVER>>(md, ls, lsolver);
+ std::unique_ptr<nonlin_model_pb<PLSOLVER>> mdpb;
+ if (md.has_internal_variables())
+ mdpb = std::make_unique<nonlin_condensed_model_pb<PLSOLVER>>(md, ls,
lsolver);
+ else
+ mdpb = std::make_unique<nonlin_model_pb<PLSOLVER>>(md, ls, lsolver);
if (dynamic_cast<newton_search_with_step_control *>(&ls))
Newton_with_step_control(*mdpb, iter);
else
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index ee9ce16..3180f5a 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -292,10 +292,15 @@ namespace getfem {
return it->second.v_num_data[niter];
}
- size_type model::nb_dof() const {
+ size_type model::nb_dof(bool with_internal) const {
context_check();
if (act_size_to_be_done) actualize_sizes();
- return complex_version ? gmm::vect_size(crhs) : gmm::vect_size(rrhs);
+ if (complex_version)
+ return gmm::vect_size(crhs);
+ else if (with_internal && gmm::vect_size(full_rrhs))
+ return gmm::vect_size(full_rrhs);
+ else
+ return gmm::vect_size(rrhs);
}
void model::resize_global_system() const {
@@ -333,6 +338,17 @@ namespace getfem {
gmm::resize(rrhs, primary_size);
}
+ if (full_size > primary_size) {
+ GMM_ASSERT1(has_internal_variables(), "Internal error");
+ gmm::resize(internal_rTM, full_size-primary_size, primary_size);
+ gmm::resize(full_rrhs, full_size);
+ gmm::resize(internal_sol, full_size-primary_size);
+ } else {
+ GMM_ASSERT1(!(has_internal_variables()), "Internal error");
+ gmm::resize(internal_rTM, 0, 0);
+ full_rrhs.clear();
+ }
+
for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib)
for (const term_description &term : bricks[ib].tlist)
if (term.is_global) {
@@ -690,8 +706,10 @@ namespace getfem {
bool firstvar(true);
for (const auto &v : variables) {
if (v.second.is_variable) {
+ const model_real_plain_vector &rhs = v.second.is_internal
+ ? full_rrhs : rrhs;
const gmm::sub_interval &II = interval_of_variable(v.first);
- scalar_type res = gmm::vect_norm2(gmm::sub_vector(rrhs, II));
+ scalar_type res = gmm::vect_norm2(gmm::sub_vector(rhs, II));
if (!firstvar) cout << ", ";
ost << "res_" << v.first << "= " << std::setw(11) << res;
firstvar = false;
@@ -794,6 +812,12 @@ namespace getfem {
act_size_to_be_done = true;
}
+ void model::add_internal_im_variable(const std::string &name,
+ const im_data &imd) {
+ add_im_variable(name, imd);
+ variables[name].is_internal = true;
+ }
+
void model::add_im_data(const std::string &name, const im_data &imd,
size_type niter) {
check_name_validity(name);
@@ -2301,11 +2325,14 @@ namespace getfem {
-
-
-
void model::assembly(build_version version) {
+ GMM_ASSERT1(version != BUILD_ON_DATA_CHANGE,
+ "Invalid assembly version BUILD_ON_DATA_CHANGE");
+ GMM_ASSERT1(version != BUILD_WITH_LIN,
+ "Invalid assembly version BUILD_WITH_LIN");
+ GMM_ASSERT1(version != BUILD_WITH_INTERNAL,
+ "Invalid assembly version BUILD_WITH_INTERNAL");
#if GETFEM_PARA_LEVEL > 1
double t_ref = MPI_Wtime();
#endif
@@ -2598,54 +2625,115 @@ namespace getfem {
if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
}
+ if (gmm::vect_size(full_rrhs)) { // i.e. if has_internal_variables()
+ gmm::sub_interval IP(0,gmm::vect_size(rrhs));
+ gmm::copy(rrhs, gmm::sub_vector(full_rrhs, IP));
+ }
+
// Generic expressions
if (generic_expressions.size()) {
- model_real_plain_vector residual;
- if (version & BUILD_RHS) gmm::resize(residual, gmm::vect_size(rrhs));
-
- { //need parentheses for constructor/destructor semantics of distro
- accumulated_distro<decltype(rrhs)> residual_distributed(residual);
- accumulated_distro<decltype(rTM)> tangent_matrix_distributed(rTM);
+ GMM_ASSERT1(!is_complex(), "to be done");
- /*running the assembly in parallel*/
- GETFEM_OMP_PARALLEL(
- if (version & BUILD_RHS)
- GMM_TRACE2("Global generic assembly RHS");
- if (version & BUILD_MATRIX)
- GMM_TRACE2("Global generic assembly tangent term");
+ if (version & BUILD_RHS)
+ GMM_TRACE2("Global generic assembly RHS");
+ if (version & BUILD_MATRIX)
+ GMM_TRACE2("Global generic assembly tangent term");
+ // auxilliary lambda function
+ auto add_assignments_and_expressions_to_workspace =
+ [&](ga_workspace &workspace)
+ {
+ for (const auto &ad : assignments)
+ workspace.add_assignment_expression
+ (ad.varname, ad.expr, ad.region, ad.order, ad.before);
+ for (const auto &ge : generic_expressions)
+ workspace.add_expression(ge.expr, ge.mim, ge.region,
+ 2, ge.secondary_domain);
+ };
+
+ const bool with_internal = version & BUILD_WITH_INTERNAL
+ && has_internal_variables();
+ model_real_sparse_matrix intern_mat; // temp for extracting condensation
info
+ model_real_plain_vector res0, // holds the original RHS
+ res1; // holds the condensed RHS
+
+ size_type full_size = gmm::vect_size(full_rrhs),
+ primary_size = gmm::vect_size(rrhs);
+
+ if ((version & BUILD_RHS) || (version & BUILD_MATRIX && with_internal))
+ gmm::resize(res0, with_internal ? full_size : primary_size);
+ if (version & BUILD_MATRIX && with_internal)
+ gmm::resize(res1, full_size);
+
+ if (version & BUILD_MATRIX) {
+ if (with_internal) {
+ gmm::resize(intern_mat, full_size, primary_size);
+ gmm::resize(res1, full_size);
+ }
+ accumulated_distro<decltype(rTM)> tangent_matrix_distro(rTM);
+ accumulated_distro<decltype(intern_mat)> intern_mat_distro(intern_mat);
+ accumulated_distro<model_real_plain_vector> res1_distro(res1);
+
+ if (version & BUILD_RHS) { // both BUILD_RHS & BUILD_MATRIX
+ gmm::resize(res0, with_internal ? full_size : primary_size);
+ accumulated_distro<model_real_plain_vector> res0_distro(res0);
+ GETFEM_OMP_PARALLEL( // running the assembly in parallel
ga_workspace workspace(*this);
-
- for (const auto &ad : assignments)
- workspace.add_assignment_expression
- (ad.varname, ad.expr, ad.region, ad.order, ad.before);
-
- for (const auto &ge : generic_expressions)
- workspace.add_expression(ge.expr, ge.mim, ge.region,
- 2, ge.secondary_domain);
-
- if (version & BUILD_RHS) {
- if (is_complex()) {
- GMM_ASSERT1(false, "to be done");
- } else {
- workspace.set_assembled_vector(residual_distributed);
- workspace.assembly(1);
- }
+ add_assignments_and_expressions_to_workspace(workspace);
+ workspace.set_assembled_vector(res0_distro);
+ workspace.assembly(1, with_internal);
+ if (with_internal) { // Condensation reads from/writes to rhs
+ gmm::copy(res0_distro.get(), res1_distro.get());
+ gmm::add(gmm::scaled(full_rrhs, scalar_type(-1)),
+ res1_distro.get()); // TIC: initial value residual=-rhs
+ workspace.set_assembled_vector(res1_distro);
+ workspace.set_internal_coupling_matrix(intern_mat_distro);
}
-
- if (version & BUILD_MATRIX) {
- if (is_complex()) {
- GMM_ASSERT1(false, "to be done");
- } else {
- workspace.set_assembled_matrix(tangent_matrix_distributed);
- workspace.assembly(2);
- }
+ workspace.set_assembled_matrix(tangent_matrix_distro);
+ workspace.assembly(2, with_internal);
+ if (with_internal) // TOC: diff from initial value, hack to make
OMP work
+ gmm::add(full_rrhs, res1_distro.get());
+ ) // end GETFEM_OMP_PARALLEL
+ } // end of res0_distro scope
+ else { // only BUILD_MATRIX
+ GETFEM_OMP_PARALLEL( // running the assembly in parallel
+ ga_workspace workspace(*this);
+ add_assignments_and_expressions_to_workspace(workspace);
+ if (with_internal) { // Condensation reads from/writes to rhs
+ gmm::copy(gmm::scaled(full_rrhs, scalar_type(-1)),
+ res1_distro.get()); // TIC: initial value residual=-rhs
+ workspace.set_assembled_vector(res1_distro);
+ workspace.set_internal_coupling_matrix(intern_mat_distro);
}
- )
- } //end of distro scope
+ workspace.set_assembled_matrix(tangent_matrix_distro);
+ workspace.assembly(2, with_internal);
+ if (with_internal) // TOC: diff from initial value, hack to make
OMP work
+ gmm::add(full_rrhs, res1_distro.get());
+ ) // end GETFEM_OMP_PARALLEL
+ }
+ } // end of tangent_matrix_distro, intern_mat_distro, res1_distro scope
+ else if (version & BUILD_RHS) {
+ gmm::resize(res0, with_internal ? full_size : primary_size);
+ accumulated_distro<model_real_plain_vector> res0_distro(res0);
+ GETFEM_OMP_PARALLEL( // running the assembly in parallel
+ ga_workspace workspace(*this);
+ add_assignments_and_expressions_to_workspace(workspace);
+ workspace.set_assembled_vector(res0_distro);
+ workspace.assembly(1, with_internal);
+ ) // end GETFEM_OMP_PARALLEL
+ } // end of res0_distro scope
if (version & BUILD_RHS)
- gmm::add(gmm::scaled(residual, scalar_type(-1)), rrhs);
+ gmm::add(gmm::scaled(res0, scalar_type(-1)), rrhs);
+
+ if (version & BUILD_MATRIX && with_internal) {
+ gmm::scale(res1, scalar_type(-1)); // from residual to rhs
+ gmm::sub_interval IP(0, primary_size),
+ II(primary_size, full_size-primary_size);
+ gmm::copy(gmm::sub_matrix(intern_mat, II, IP), internal_rTM); // -->
internal_rTM
+ gmm::add(gmm::sub_vector(res1, IP), rrhs); // -->
rrhs
+ gmm::copy(gmm::sub_vector(res1, II), internal_sol); // -->
internal_sol
+ }
}
// Post simplification for dof constraints