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

Fix multiple linear/nonlinear solves #241

Merged
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
1 change: 0 additions & 1 deletion include/model_variables.h
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,6 @@ class modelResidual
struct variable_info
{
bool is_scalar;
unsigned int variable_index;
unsigned int global_var_index;
dealii::EvaluationFlags::EvaluationFlags evaluation_flags;
dealii::EvaluationFlags::EvaluationFlags residual_flags;
Expand Down
46 changes: 25 additions & 21 deletions include/variableContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,14 @@
#ifndef VARIBLECONTAINER_H
#define VARIBLECONTAINER_H

#include <deal.II/base/exceptions.h>
#include <deal.II/lac/vector.h>
#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>

#include <boost/unordered_map.hpp>

#include "userInputParameters.h"

template <int dim, int degree, typename T>
Expand All @@ -19,15 +22,16 @@ class variableContainer

// Standard contructor, used for most situations
variableContainer(const dealii::MatrixFree<dim, double> &data,
std::vector<variable_info> _varInfoList,
std::vector<variable_info> _varChangeInfoList);
const std::vector<variable_info> &_varInfoList,
const std::vector<variable_info> &_varChangeInfoList);

variableContainer(const dealii::MatrixFree<dim, double> &data,
std::vector<variable_info> _varInfoList);
const std::vector<variable_info> &_varInfoList);
// Nonstandard constructor, used when only one index of "data" should be used,
// use with care!
variableContainer(const dealii::MatrixFree<dim, double> &data,
std::vector<variable_info> _varInfoList,
unsigned int fixed_index);
const std::vector<variable_info> &_varInfoList,
const unsigned int &fixed_index);

// Methods to get the value/grad/hess in the residual method (this is how the
// user gets these values in equations.h)
Expand Down Expand Up @@ -90,11 +94,6 @@ class variableContainer
reinit_and_eval_change_in_solution(const vectorType &src,
unsigned int cell,
unsigned int var_being_solved);
void
reinit_and_eval_LHS(const vectorType &src,
const std::vector<vectorType *> solutionSet,
unsigned int cell,
unsigned int var_being_solved);

// Only initialize the FEEvaluation object for each variable (used for
// post-processing)
Expand All @@ -111,29 +110,34 @@ class variableContainer
// The quadrature point index, a method to get the number of quadrature points
// per cell, and a method to get the xyz coordinates for the quadrature point
unsigned int q_point;

unsigned int
get_num_q_points();
get_num_q_points() const;

dealii::Point<dim, T>
get_q_point_location();
get_q_point_location() const;

private:
// The number of variables
unsigned int num_var;

// Vectors of the actual FEEvaluation objects for each active variable, split
// into scalar variables and vector variables for type reasons
std::vector<dealii::FEEvaluation<dim, degree, degree + 1, 1, double>> scalar_vars;
std::vector<dealii::FEEvaluation<dim, degree, degree + 1, dim, double>> vector_vars;
using scalar_FEEval = dealii::FEEvaluation<dim, degree, degree + 1, 1, double>;
using vector_FEEval = dealii::FEEvaluation<dim, degree, degree + 1, dim, double>;

boost::unordered_map<unsigned int, std::unique_ptr<scalar_FEEval>> scalar_vars_map;
boost::unordered_map<unsigned int, std::unique_ptr<vector_FEEval>> vector_vars_map;

std::vector<dealii::FEEvaluation<dim, degree, degree + 1, 1, double>>
scalar_change_in_vars;
std::vector<dealii::FEEvaluation<dim, degree, degree + 1, dim, double>>
vector_change_in_vars;
boost::unordered_map<unsigned int, std::unique_ptr<scalar_FEEval>>
scalar_change_in_vars_map;
boost::unordered_map<unsigned int, std::unique_ptr<vector_FEEval>>
vector_change_in_vars_map;

// Object containing some information about each variable (indices, whether
// the val/grad/hess is needed, etc)
std::vector<variable_info> varInfoList;
std::vector<variable_info> varChangeInfoList;

// The number of variables
unsigned int num_var;
};

#endif
1 change: 0 additions & 1 deletion src/matrixfree/computeLHS.cc
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@ MatrixFreePDE<dim, degree>::getLHS(
for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
{
// Initialize, read DOFs, and set evaulation flags for each variable
// variable_list.reinit_and_eval_LHS(src,solutionSet,cell,currentFieldIndex);
variable_list.reinit_and_eval(solutionSet, cell);
variable_list.reinit_and_eval_change_in_solution(src, cell, currentFieldIndex);

Expand Down
142 changes: 13 additions & 129 deletions src/userInputParameters/loadVariableAttributes.cc
Original file line number Diff line number Diff line change
Expand Up @@ -50,8 +50,6 @@ userInputParameters<dim>::loadVariableAttributes(
}
}
varInfoListExplicitRHS.reserve(num_var_explicit_RHS);
unsigned int scalar_var_index = 0;
unsigned int vector_var_index = 0;
for (unsigned int i = 0; i < number_of_variables; i++)
{
variable_info varInfo;
Expand All @@ -64,28 +62,9 @@ userInputParameters<dim>::loadVariableAttributes(

varInfo.global_var_index = i;

!(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
? varInfo.var_needed = true
: varInfo.var_needed = false;
varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);

if (var_type[i] == SCALAR)
{
varInfo.is_scalar = true;
if (varInfo.var_needed)
{
varInfo.variable_index = scalar_var_index;
scalar_var_index++;
}
}
else
{
varInfo.is_scalar = false;
if (varInfo.var_needed)
{
varInfo.variable_index = vector_var_index;
vector_var_index++;
}
}
varInfo.is_scalar = var_type[i] == SCALAR;

varInfoListExplicitRHS.push_back(varInfo);
}
Expand All @@ -101,8 +80,6 @@ userInputParameters<dim>::loadVariableAttributes(
}
}
varInfoListNonexplicitRHS.reserve(num_var_nonexplicit_RHS);
scalar_var_index = 0;
vector_var_index = 0;
for (unsigned int i = 0; i < number_of_variables; i++)
{
variable_info varInfo;
Expand All @@ -115,28 +92,9 @@ userInputParameters<dim>::loadVariableAttributes(

varInfo.global_var_index = i;

!(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
? varInfo.var_needed = true
: varInfo.var_needed = false;
varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);

if (var_type[i] == SCALAR)
{
varInfo.is_scalar = true;
if (varInfo.var_needed)
{
varInfo.variable_index = scalar_var_index;
scalar_var_index++;
}
}
else
{
varInfo.is_scalar = false;
if (varInfo.var_needed)
{
varInfo.variable_index = vector_var_index;
vector_var_index++;
}
}
varInfo.is_scalar = var_type[i] == SCALAR;

varInfoListNonexplicitRHS.push_back(varInfo);
}
Expand All @@ -153,8 +111,6 @@ userInputParameters<dim>::loadVariableAttributes(
}

varInfoListLHS.reserve(num_var_LHS);
scalar_var_index = 0;
vector_var_index = 0;
for (unsigned int i = 0; i < number_of_variables; i++)
{
variable_info varInfo;
Expand All @@ -167,35 +123,14 @@ userInputParameters<dim>::loadVariableAttributes(

varInfo.global_var_index = i;

!(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
? varInfo.var_needed = true
: varInfo.var_needed = false;
varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);

if (var_type[i] == SCALAR)
{
varInfo.is_scalar = true;
if (varInfo.var_needed)
{
varInfo.variable_index = scalar_var_index;
scalar_var_index++;
}
}
else
{
varInfo.is_scalar = false;
if (varInfo.var_needed)
{
varInfo.variable_index = vector_var_index;
vector_var_index++;
}
}
varInfo.is_scalar = var_type[i] == SCALAR;

varInfoListLHS.push_back(varInfo);
}

varChangeInfoListLHS.reserve(num_var_LHS);
scalar_var_index = 0;
vector_var_index = 0;
for (unsigned int i = 0; i < number_of_variables; i++)
{
variable_info varInfo;
Expand All @@ -209,37 +144,16 @@ userInputParameters<dim>::loadVariableAttributes(

varInfo.global_var_index = i;

!(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
? varInfo.var_needed = true
: varInfo.var_needed = false;
varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);

if (var_type[i] == SCALAR)
{
varInfo.is_scalar = true;
if (varInfo.var_needed)
{
varInfo.variable_index = scalar_var_index;
scalar_var_index++;
}
}
else
{
varInfo.is_scalar = false;
if (varInfo.var_needed)
{
varInfo.variable_index = vector_var_index;
vector_var_index++;
}
}
varInfo.is_scalar = var_type[i] == SCALAR;

varChangeInfoListLHS.push_back(varInfo);
}

// Load variable information for postprocessing
// First, the info list for the base field variables
pp_baseVarInfoList.reserve(number_of_variables);
scalar_var_index = 0;
vector_var_index = 0;
for (unsigned int i = 0; i < number_of_variables; i++)
{
variable_info varInfo;
Expand All @@ -249,28 +163,9 @@ userInputParameters<dim>::loadVariableAttributes(

varInfo.global_var_index = i;

!(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing)
? varInfo.var_needed = true
: varInfo.var_needed = false;
varInfo.var_needed = !(varInfo.evaluation_flags & dealii::EvaluationFlags::nothing);

if (var_type[i] == SCALAR)
{
varInfo.is_scalar = true;
if (varInfo.var_needed)
{
varInfo.variable_index = scalar_var_index;
scalar_var_index++;
}
}
else
{
varInfo.is_scalar = false;
if (varInfo.var_needed)
{
varInfo.variable_index = vector_var_index;
vector_var_index++;
}
}
varInfo.is_scalar = var_type[i] == SCALAR;

pp_baseVarInfoList.push_back(varInfo);
}
Expand All @@ -292,8 +187,6 @@ userInputParameters<dim>::loadVariableAttributes(

// The info list for the postprocessing field variables
pp_varInfoList.reserve(pp_number_of_variables);
scalar_var_index = 0;
vector_var_index = 0;
for (unsigned int i = 0; i < pp_number_of_variables; i++)
{
variable_info varInfo;
Expand All @@ -303,18 +196,9 @@ userInputParameters<dim>::loadVariableAttributes(
variable_attributes.equation_dependency_parser.eval_flags_residual_postprocess[i];

varInfo.global_var_index = i;
if (pp_var_type[i] == SCALAR)
{
varInfo.is_scalar = true;
varInfo.variable_index = scalar_var_index;
scalar_var_index++;
}
else
{
varInfo.is_scalar = false;
varInfo.variable_index = vector_var_index;
vector_var_index++;
}

varInfo.is_scalar = pp_var_type[i] == SCALAR;

pp_varInfoList.push_back(varInfo);
}
}
Expand Down
Loading