Skip to content

[WIP] Allow specifying element range in project_vector #4263

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.

Sign up for GitHub

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

Draft
wants to merge 1 commit into
base: devel
Choose a base branch
from
Draft

[WIP] Allow specifying element range in project_vector #4263

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 15 additions & 7 deletions include/systems/system.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@
#include <string_view>
#include <vector>
#include <memory>
#include <optional>

// This define may be useful in --disable-optional builds when it is
// possible that libmesh will not have any solvers available.
Expand Down Expand Up @@ -561,7 +562,8 @@ class System : public ReferenceCountedObject<System>,
void project_vector (NumericVector<Number> & new_vector,
FunctionBase<Number> * f,
FunctionBase<Gradient> * g = nullptr,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/**
* Projects arbitrary functions onto a vector of degree of freedom
Expand All @@ -577,7 +579,8 @@ class System : public ReferenceCountedObject<System>,
void project_vector (NumericVector<Number> & new_vector,
FEMFunctionBase<Number> * f,
FEMFunctionBase<Gradient> * g = nullptr,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/**
* Projects arbitrary functions onto a vector of degree of freedom
Expand All @@ -594,7 +597,8 @@ class System : public ReferenceCountedObject<System>,
GradientFunctionPointer gptr,
const Parameters & parameters,
NumericVector<Number> & new_vector,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/**
* Projects arbitrary boundary functions onto a vector of degree of
Expand Down Expand Up @@ -651,7 +655,8 @@ class System : public ReferenceCountedObject<System>,
NumericVector<Number> & new_vector,
FunctionBase<Number> * f,
FunctionBase<Gradient> * g = nullptr,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/**
* Projects arbitrary boundary functions onto a vector of degree of
Expand All @@ -674,7 +679,8 @@ class System : public ReferenceCountedObject<System>,
GradientFunctionPointer gptr,
const Parameters & parameters,
NumericVector<Number> & new_vector,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/**
* \returns The system number.
Expand Down Expand Up @@ -2005,7 +2011,8 @@ class System : public ReferenceCountedObject<System>,
* primal constraints if is_adjoint is non-negative.
*/
void project_vector (NumericVector<Number> &,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/**
* Projects the vector defined on the old mesh onto the
Expand All @@ -2017,7 +2024,8 @@ class System : public ReferenceCountedObject<System>,
*/
void project_vector (const NumericVector<Number> &,
NumericVector<Number> &,
int is_adjoint = -1) const;
int is_adjoint = -1,
std::optional<ConstElemRange> active_local_range = std::nullopt) const;

/*
* If we have e.g. a element space constrained by spline values, we
Expand Down
73 changes: 46 additions & 27 deletions src/systems/system_projection.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -245,15 +245,16 @@ public:
// ------------------------------------------------------------
// System implementation
void System::project_vector (NumericVector<Number> & vector,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
// Create a copy of the vector, which currently
// contains the old data.
std::unique_ptr<NumericVector<Number>>
old_vector (vector.clone());

// Project the old vector to the new vector
this->project_vector (*old_vector, vector, is_adjoint);
this->project_vector (*old_vector, vector, is_adjoint, active_local_range);
}


Expand All @@ -264,7 +265,8 @@ void System::project_vector (NumericVector<Number> & vector,
*/
void System::project_vector (const NumericVector<Number> & old_v,
NumericVector<Number> & new_v,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
LOG_SCOPE ("project_vector(old,new)", "System");

Expand All @@ -285,9 +287,12 @@ void System::project_vector (const NumericVector<Number> & old_v,
std::unique_ptr<NumericVector<Number>> local_old_vector_built;
const NumericVector<Number> * old_vector_ptr = nullptr;

ConstElemRange active_local_elem_range
(this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end());
if (!active_local_range)
{
active_local_range.emplace
(this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end());
}

// If the old vector was uniprocessor, make the new
// vector uniprocessor
Expand All @@ -304,7 +309,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
{
// Build a send list for efficient localization
BuildProjectionList projection_list(*this);
Threads::parallel_reduce (active_local_elem_range,
Threads::parallel_reduce (active_local_range.value(),
projection_list);

// Create a sorted, unique send_list
Expand All @@ -328,7 +333,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
{
// Build a send list for efficient localization
BuildProjectionList projection_list(*this);
Threads::parallel_reduce (active_local_elem_range,
Threads::parallel_reduce (active_local_range.value(),
projection_list);

// Create a sorted, unique send_list
Expand Down Expand Up @@ -389,7 +394,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
g(*this, old_vector, &regular_vars);

FEMProjector projector(*this, f, &g, setter, regular_vars);
projector.project(active_local_elem_range);
projector.project(active_local_range.value());
}

if (!vector_vars.empty())
Expand All @@ -403,7 +408,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
OldSolutionValue<Tensor, &FEMContext::point_gradient> g_vector(*this, old_vector, &vector_vars);

FEMVectorProjector vector_projector(*this, f_vector, &g_vector, setter, vector_vars);
vector_projector.project(active_local_elem_range);
vector_projector.project(active_local_range.value());
}

// Copy the SCALAR dofs from old_vector to new_vector
Expand Down Expand Up @@ -1068,11 +1073,12 @@ void System::project_vector (ValueFunctionPointer fptr,
GradientFunctionPointer gptr,
const Parameters & function_parameters,
NumericVector<Number> & new_vector,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
WrappedFunction<Number> f(*this, fptr, &function_parameters);
WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
this->project_vector(new_vector, &f, &g, is_adjoint);
this->project_vector(new_vector, &f, &g, is_adjoint, active_local_range);
}

/**
Expand All @@ -1082,7 +1088,8 @@ void System::project_vector (ValueFunctionPointer fptr,
void System::project_vector (NumericVector<Number> & new_vector,
FunctionBase<Number> * f,
FunctionBase<Gradient> * g,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
LOG_SCOPE ("project_vector(FunctionBase)", "System");

Expand All @@ -1094,10 +1101,10 @@ void System::project_vector (NumericVector<Number> & new_vector,
{
WrappedFunctor<Gradient> g_fem(*g);

this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint);
this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint, active_local_range);
}
else
this->project_vector(new_vector, &f_fem, nullptr, is_adjoint);
this->project_vector(new_vector, &f_fem, nullptr, is_adjoint, active_local_range);
}


Expand All @@ -1108,15 +1115,19 @@ void System::project_vector (NumericVector<Number> & new_vector,
void System::project_vector (NumericVector<Number> & new_vector,
FEMFunctionBase<Number> * f,
FEMFunctionBase<Gradient> * g,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
LOG_SCOPE ("project_fem_vector()", "System");

libmesh_assert (f);

ConstElemRange active_local_range
(this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end() );
if (!active_local_range)
{
active_local_range.emplace
(this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end());
}

VectorSetAction<Number> setter(new_vector);

Expand All @@ -1137,12 +1148,12 @@ void System::project_vector (NumericVector<Number> & new_vector,
FEMFunctionWrapper<Gradient> gw(*g);

FEMProjector projector(*this, fw, &gw, setter, vars);
projector.project(active_local_range);
projector.project(active_local_range.value());
}
else
{
FEMProjector projector(*this, fw, nullptr, setter, vars);
projector.project(active_local_range);
projector.project(active_local_range.value());
}

// Also, load values into the SCALAR dofs
Expand Down Expand Up @@ -1197,7 +1208,7 @@ void System::project_vector (NumericVector<Number> & new_vector,
{
// Look for a spline node: a NodeElem with a rational variable
// on it.
for (auto & elem : active_local_range)
for (auto & elem : active_local_range.value())
if (elem->type() == NODEELEM)
for (auto rational_var : rational_vars)
if (rational_var->active_on_subdomain(elem->subdomain_id()))
Expand Down Expand Up @@ -1274,12 +1285,13 @@ void System::boundary_project_vector (const std::set<boundary_id_type> & b,
GradientFunctionPointer gptr,
const Parameters & function_parameters,
NumericVector<Number> & new_vector,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
WrappedFunction<Number> f(*this, fptr, &function_parameters);
WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
this->boundary_project_vector(b, variables, new_vector, &f, &g,
is_adjoint);
is_adjoint, active_local_range);
}

/**
Expand All @@ -1291,13 +1303,20 @@ void System::boundary_project_vector (const std::set<boundary_id_type> & b,
NumericVector<Number> & new_vector,
FunctionBase<Number> * f,
FunctionBase<Gradient> * g,
int is_adjoint) const
int is_adjoint,
std::optional<ConstElemRange> active_local_range) const
{
LOG_SCOPE ("boundary_project_vector()", "System");

if (!active_local_range)
{
active_local_range.emplace
(this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end());
}

Threads::parallel_for
(ConstElemRange (this->get_mesh().active_local_elements_begin(),
this->get_mesh().active_local_elements_end() ),
(active_local_range.value(),
BoundaryProjectSolution(b, variables, *this, f, g,
this->get_equation_systems().parameters,
new_vector)
Expand Down