Skip to content

Commit

Permalink
Add temperature extrapolation (#994)
Browse files Browse the repository at this point in the history 
Description
This PR adds the extrapolation in time of temperature and temperature gradient value for source terms in the Navier-Stokes equations requiring them (Marangoni and evaporation recoil force terms).
How Has This Been Tested?
All tests have passed
Comment
Velocity extrapolation has been moved to include/solvers/vof_scratch_data.h. Solution extrapolations are now located in ScratchData to avoid constraint on assemblers order.
  • Loading branch information
@AmishgaAlphonius
AmishgaAlphonius authored Jan 24, 2024
1 parent c39dbed commit 6dc4461
Show file tree
Hide file tree
Showing 8 changed files with 152 additions and 70 deletions.
102 changes: 74 additions & 28 deletions include/solvers/navier_stokes_scratch_data.h
View file
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
View file
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
View file
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
View file
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,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 @@ -211,6 +211,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 @@ -340,6 +341,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 @@ -431,9 +433,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

0 comments on commit 6dc4461

Please sign in to comment.