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

Add temperature extrapolation #994

Merged
merged 2 commits into from
Jan 24, 2024
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
102 changes: 74 additions & 28 deletions include/solvers/navier_stokes_scratch_data.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <core/parameters.h>
#include <core/physical_property_model.h>
#include <core/rheological_model.h>
#include <core/time_integration_utilities.h>

#include <solvers/physical_properties_manager.h>
#include <solvers/vof_filter.h>
Expand Down Expand Up @@ -76,19 +77,31 @@ class NavierStokesScratchData
* necessary memory for all member variables. However, it does not do any
* evaluation, since this needs to be done at the cell level.
*
* @param simulation_control The SimulationControl object that holds
* information related to the control of the steady-state or transient
* simulation. This is used to extrapolate auxiliary physics solutions in time
* for transient simulation.
*
* @param properties_manager The PhysicalPropertiesManager object that stores
* physical property models.
*
* @param fe The FESystem used to solve the Navier-Stokes equations
*
* @param quadrature The quadrature to use for the assembly
*
* @param mapping The mapping of the domain in which the Navier-Stokes equations are solved
* @param mapping The mapping of the domain in which the Navier-Stokes
* equations are solved
*
*/
NavierStokesScratchData(PhysicalPropertiesManager &properties_manager,
const FESystem<dim> &fe,
const Quadrature<dim> &quadrature,
const Mapping<dim> &mapping,
const Quadrature<dim - 1> &face_quadrature)
: properties_manager(properties_manager)
NavierStokesScratchData(
const std::shared_ptr<SimulationControl> &simulation_control,
const PhysicalPropertiesManager &properties_manager,
const FESystem<dim> &fe,
const Quadrature<dim> &quadrature,
const Mapping<dim> &mapping,
const Quadrature<dim - 1> &face_quadrature)
: simulation_control(simulation_control)
, properties_manager(properties_manager)
, fe_values(mapping,
fe,
quadrature,
Expand Down Expand Up @@ -117,19 +130,16 @@ class NavierStokesScratchData

/**
* @brief Copy Constructor. Same as the main constructor.
* This constructor only uses the other scratch to build the FeValues, it
* This constructor only uses the other scratch to build the FeValues, it
* does not copy the content of the other scratch into itself since, by
* definition of the WorkStream mechanism it is assumed that the content of
* the scratch will be reset on a cell basis.
*
* @param fe The FESystem used to solve the Navier-Stokes equations
*
* @param quadrature The quadrature to use for the assembly
*
* @param mapping The mapping of the domain in which the Navier-Stokes equations are solved
* @param sd The scratch data to be copied
*/
NavierStokesScratchData(const NavierStokesScratchData<dim> &sd)
: properties_manager(sd.properties_manager)
: simulation_control(sd.simulation_control)
, properties_manager(sd.properties_manager)
, fe_values(sd.fe_values.get_mapping(),
sd.fe_values.get_fe(),
sd.fe_values.get_quadrature(),
Expand Down Expand Up @@ -578,7 +588,8 @@ class NavierStokesScratchData
}

/**
* @brief enable_void_fraction Enables the collection of the void fraction data by the scratch
* @brief enable_void_fraction Enables the collection of the void fraction
* data by the scratch
*
* @param fe FiniteElement associated with the void fraction
*
Expand Down Expand Up @@ -642,15 +653,17 @@ class NavierStokesScratchData
const unsigned int n_global_max_particles_per_cell,
const bool enable_void_fraction_interpolation);

/** @brief Calculate the content of the scratch for the particle fluid interactions
/** @brief Calculates the content of the scratch for the particle fluid
* interactions
*
* @param cell The cell over which the assembly is being carried.
* This cell must be compatible with the void fraction FE and not the
* Navier-Stokes FE
*
* @param current_solution The present value of the solution for [epsilon]
*
* @param previous_solutions The solutions at the previous time steps for [epsilon]
* @param previous_solutions The solutions at the previous time steps for
* [epsilon]
*
*/

Expand Down Expand Up @@ -808,7 +821,7 @@ class NavierStokesScratchData
// If particles are in the cell, gather the rest of the
// information

// Create local vector that will be use to spawn an in-situ quadrature to
// Create local vector that will be used to spawn an in-situ quadrature to
// interpolate at the location of the particles
std::vector<Point<dim>> particle_reference_location(number_of_particles);
std::vector<double> particle_weights(number_of_particles, 1);
Expand Down Expand Up @@ -878,23 +891,26 @@ class NavierStokesScratchData


/** @brief Reinitialize the content of the scratch for the heat transfer
* auxiliary physic
*
* @param cell The cell over which the assembly is being carried.
* This cell must be compatible with the heat transfer FE and not the
* Navier-Stokes FE
*
* @param current_solution The present value of the solution for temperature
*
* @param previous_solutions The solutions at the previous time steps for
* temperature
* @param previous_solutions Vector of \f$n\f$ @p VectorType containers of
* previous temperature solutions. \f$n\f$ depends on the BDF scheme selected
* for time-stepping.
*
*/

template <typename VectorType>
void
reinit_heat_transfer(
const typename DoFHandler<dim>::active_cell_iterator &cell,
const VectorType &current_solution)
const VectorType &current_solution,
const std::vector<VectorType> &previous_solutions)
{
this->fe_values_temperature->reinit(cell);

Expand All @@ -905,6 +921,32 @@ class NavierStokesScratchData
// Gather temperature gradient
this->fe_values_temperature->get_function_gradients(
current_solution, this->temperature_gradients);

// Extrapolate temperature and temperature gradient to t+dt using the BDF
// if the simulation is transient
const auto method = this->simulation_control->get_assembly_method();
if (is_bdf(method))
{
// Gather previous temperature and temperature gradient values
for (unsigned int p = 0; p < previous_solutions.size(); ++p)
{
fe_values_temperature->get_function_values(
previous_solutions[p], this->previous_temperature_values[p]);
fe_values_temperature->get_function_gradients(
previous_solutions[p], this->previous_temperature_gradients[p]);
}
// Extrapolate temperature and temperature gradient
std::vector<double> time_vector =
this->simulation_control->get_simulation_times();
bdf_extrapolate(time_vector,
previous_temperature_values,
number_of_previous_solutions(method),
temperature_values);
bdf_extrapolate(time_vector,
previous_temperature_gradients,
number_of_previous_solutions(method),
temperature_gradients);
}
}

/**
Expand Down Expand Up @@ -973,8 +1015,11 @@ class NavierStokesScratchData
void
calculate_physical_properties();

// For auxiliary physics solution extrapolation
const std::shared_ptr<SimulationControl> simulation_control;

// Physical properties
PhysicalPropertiesManager properties_manager;
const PhysicalPropertiesManager properties_manager;
std::map<field, std::vector<double>> fields;
std::vector<double> density;
double density_ref;
Expand Down Expand Up @@ -1049,7 +1094,6 @@ class NavierStokesScratchData
std::vector<std::vector<double>> phi_p;
std::vector<std::vector<Tensor<1, dim>>> grad_phi_p;


/**
* Scratch component for the VOF auxiliary physics
*/
Expand Down Expand Up @@ -1110,10 +1154,12 @@ class NavierStokesScratchData
/**
* Scratch component for the heat transfer
*/
bool gather_temperature;
unsigned int n_dofs_heat_transfer;
std::vector<double> temperature_values;
std::vector<Tensor<1, dim>> temperature_gradients;
bool gather_temperature;
unsigned int n_dofs_heat_transfer;
std::vector<double> temperature_values;
std::vector<std::vector<double>> previous_temperature_values;
std::vector<Tensor<1, dim>> temperature_gradients;
std::vector<std::vector<Tensor<1, dim>>> previous_temperature_gradients;
// This is stored as a shared_ptr because it is only instantiated when needed
std::shared_ptr<FEValues<dim>> fe_values_temperature;

Expand All @@ -1139,7 +1185,7 @@ class NavierStokesScratchData
*/
bool is_boundary_cell;

// If a rheological model is being used for a non Newtonian flow
// If a rheological model is being used for a non-Newtonian flow
bool gather_hessian;

FEFaceValues<dim> fe_face_values;
Expand Down
57 changes: 45 additions & 12 deletions include/solvers/vof_scratch_data.h
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
*/

#include <core/multiphysics.h>
#include <core/time_integration_utilities.h>

#include <solvers/multiphysics_interface.h>

Expand Down Expand Up @@ -63,23 +64,29 @@ class VOFScratchData
* necessary memory for all member variables. However, it does not do any
* evaluation, since this needs to be done at the cell level.
*
* @param properties_manager The physical properties Manager (see physical_properties_manager.h)
* @param simulation_control The SimulationControl object that holds
* information related to the control of the steady-state or transient
* simulation. This is used to extrapolate velocity solutions in time
* for transient simulation.
*
* @param fe_vof The FESystem used to solve the VOF equations
*
* @param quadrature The quadrature to use for the assembly
*
* @param mapping The mapping of the domain in which the Navier-Stokes equations are solved
* @param mapping The mapping of the domain in which the Navier-Stokes
* equations are solved
*
* @param fe_fd The FESystem used to solve the Fluid Dynamics equations
*
*/
VOFScratchData(const PhysicalPropertiesManager properties_manager,
const FiniteElement<dim> &fe_vof,
const Quadrature<dim> &quadrature,
const Mapping<dim> &mapping,
const FiniteElement<dim> &fe_fd)
: properties_manager(properties_manager)
VOFScratchData(const std::shared_ptr<SimulationControl> &simulation_control,
const PhysicalPropertiesManager &properties_manager,
const FiniteElement<dim> &fe_vof,
const Quadrature<dim> &quadrature,
const Mapping<dim> &mapping,
const FiniteElement<dim> &fe_fd)
: simulation_control(simulation_control)
, properties_manager(properties_manager)
, fe_values_vof(mapping,
fe_vof,
quadrature,
Expand All @@ -98,10 +105,12 @@ class VOFScratchData
* definition of the WorkStream mechanism it is assumed that the content of
* the scratch will be reset on a cell basis.
*
* @param sd The scratch data
* @param sd The scratch data to be copied
*/
VOFScratchData(const VOFScratchData<dim> &sd)
: fe_values_vof(sd.fe_values_vof.get_mapping(),
: simulation_control(sd.simulation_control)
, properties_manager(sd.properties_manager)
, fe_values_vof(sd.fe_values_vof.get_mapping(),
sd.fe_values_vof.get_fe(),
sd.fe_values_vof.get_quadrature(),
update_values | update_gradients |
Expand Down Expand Up @@ -196,7 +205,13 @@ class VOFScratchData
* @param cell The cell for which the velocity is reinitialized
* This cell must be compatible with the Fluid Dynamics FE
*
* @param current_solution The present value of the solution for [u,p]
* @param current_solution The present value of the solution for \f$[u,p]\f$
*
* @param previous_solutions Vector of \f$n\f$ @p VectorType containers of
* previous fluid dynamic solutions (\f$[u,p]\f$). \f$n\f$ depends on the BDF
* scheme selected for time-stepping.
*
* @param ale ALE parameters
*
*/

Expand All @@ -220,6 +235,7 @@ class VOFScratchData
trace(this->velocity_gradient_values[q]);
}

// Gather previous velocity values
for (unsigned int p = 0; p < previous_solutions.size(); ++p)
{
fe_values_fd[velocities_fd].get_function_values(
Expand Down Expand Up @@ -250,10 +266,27 @@ class VOFScratchData
this->previous_velocity_values[p][q] -= velocity_ale;
}
}

// Extrapolate velocity to t+dt using the BDF scheme if the simulation is
// transient
const auto method = this->simulation_control->get_assembly_method();
if (is_bdf(method))
{
// Extrapolate velocity
std::vector<double> time_vector =
this->simulation_control->get_simulation_times();
bdf_extrapolate(time_vector,
previous_velocity_values,
number_of_previous_solutions(method),
velocity_values);
}
}

// For velocity solution extrapolation
const std::shared_ptr<SimulationControl> simulation_control;

// Physical properties
PhysicalPropertiesManager properties_manager;
const PhysicalPropertiesManager properties_manager;
std::map<field, std::vector<double>> fields;

// FEValues for the VOF problem
Expand Down
2 changes: 2 additions & 0 deletions source/fem-dem/gls_vans.cc
Original file line number Diff line number Diff line change
Expand Up @@ -1552,6 +1552,7 @@ GLSVANSSolver<dim>::assemble_system_matrix()
setup_assemblers();

auto scratch_data = NavierStokesScratchData<dim>(
this->simulation_control,
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
Expand Down Expand Up @@ -1657,6 +1658,7 @@ GLSVANSSolver<dim>::assemble_system_rhs()
setup_assemblers();

auto scratch_data = NavierStokesScratchData<dim>(
this->simulation_control,
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
Expand Down
11 changes: 7 additions & 4 deletions source/solvers/gd_navier_stokes.cc
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ GDNavierStokesSolver<dim>::setup_assemblers()
if (this->simulation_parameters.physical_properties_manager
.is_non_newtonian())
{
// Core assembler with Non newtonian viscosity
// Core assembler with non-Newtonian viscosity
this->assemblers.push_back(
std::make_shared<GDNavierStokesAssemblerNonNewtonianCore<dim>>(
this->simulation_control, gamma));
Expand Down Expand Up @@ -203,6 +203,7 @@ GDNavierStokesSolver<dim>::assemble_system_matrix()
setup_assemblers();

auto scratch_data = NavierStokesScratchData<dim>(
this->simulation_control,
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
Expand Down Expand Up @@ -332,6 +333,7 @@ GDNavierStokesSolver<dim>::assemble_system_rhs()
setup_assemblers();

auto scratch_data = NavierStokesScratchData<dim>(
this->simulation_control,
this->simulation_parameters.physical_properties_manager,
*this->fe,
*this->cell_quadrature,
Expand Down Expand Up @@ -423,9 +425,10 @@ GDNavierStokesSolver<dim>::assemble_local_system_rhs(
cell->index(),
dof_handler_ht);

scratch_data.reinit_heat_transfer(temperature_cell,
*this->multiphysics->get_solution(
PhysicsID::heat_transfer));
scratch_data.reinit_heat_transfer(
temperature_cell,
*this->multiphysics->get_solution(PhysicsID::heat_transfer),
*this->multiphysics->get_previous_solutions(PhysicsID::heat_transfer));
}

scratch_data.calculate_physical_properties();
Expand Down
Loading
Loading