Skip to content

Commit

Permalink
Advection works just fine with tracer now
Browse files Browse the repository at this point in the history
  • Loading branch information
@blaisb @oguevremont
blaisb authored and oguevremont committed Sep 24, 2024
1 parent 072c73a commit 5e15ceb
Show file tree
Hide file tree
Showing 4 changed files with 115 additions and 74 deletions.
21 changes: 16 additions & 5 deletions include/core/boundary_conditions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ namespace BoundaryConditions
partial_slip,
periodic,
pressure,
outlet,
// for heat transfer
outlet, // outlet is also used for tracers and eventually for other physics
// for heat transfer
noflux,
temperature,
convection_radiation,
Expand All @@ -58,7 +58,7 @@ namespace BoundaryConditions
cahn_hilliard_noflux,
cahn_hilliard_dirichlet_phase_order,
cahn_hilliard_angle_of_contact,
cahn_hilliard_free_angle,
cahn_hilliard_free_angle
};

/**
Expand Down Expand Up @@ -720,9 +720,9 @@ namespace BoundaryConditions
{
prm.declare_entry("type",
"dirichlet",
Patterns::Selection("dirichlet"),
Patterns::Selection("dirichlet|outlet"),
"Type of boundary condition for tracer"
"Choices are <function>.");
"Choices are <dirichlet|outlet>.");

prm.declare_entry("id",
Utilities::int_to_string(i_bc, 2),
Expand Down Expand Up @@ -798,6 +798,17 @@ namespace BoundaryConditions
prm.leave_subsection();
}

else if (op == "outlet")
{
this->type[i_bc] = BoundaryType::outlet;
}
else
{
AssertThrow(false,
ExcMessage("Unknown boundary condition type for tracers"));
}


this->id[i_bc] = prm.get_integer("id");
}

Expand Down
51 changes: 51 additions & 0 deletions include/solvers/tracer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,29 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
// outputs
if (simulation_parameters.timer.type == Parameters::Timer::Type::none)
this->computing_timer.disable_output();

/// Verify that we have specified as many tracer boundary conditions as
/// there are in the triangulation
// const std::vector<types::boundary_id> boundary_ids_in_triangulation =
// p_triangulation->get_boundary_ids();

/// Verify that the tracer boundary conditions are specified for all
/// boundary ids
// for (unsigned int i_bc = 0;
// i_bc < simulation_parameters.boundary_conditions_tracer.size;
// ++i_bc)
// {
// Assert(std::find(boundary_ids_in_triangulation.begin(),
// boundary_ids_in_triangulation.end(),
// simulation_parameters.boundary_conditions_tracer
// .id[i_bc]) != boundary_ids_in_triangulation.end(),
// ExcMessage(
// "The tracer boundary condition with id " +
// std::to_string(
// simulation_parameters.boundary_conditions_tracer.id[i_bc])
// +
// " is not present in the triangulation."));
// }
}

/**
Expand Down Expand Up @@ -412,6 +435,34 @@ class Tracer : public AuxiliaryPhysics<dim, GlobalVectorType>
void
write_tracer_statistics();


/**
* @brief Get the lethe boundary indicator for a given triangulation boundary while carrying the appropriate checks
*/
unsigned int
get_lethe_boundary_index(const types::boundary_id &triangulation_boundary_id)
{
// Identify which boundary condition corresponds to the boundary id. If
// this boundary condition is not identified, then exit the simulation
// instead of assuming an outlet.
const auto lethe_boundary_id = std::find(
this->simulation_parameters.boundary_conditions_tracer.id.begin(),
this->simulation_parameters.boundary_conditions_tracer.id.end(),
triangulation_boundary_id);

AssertThrow(
lethe_boundary_id !=
this->simulation_parameters.boundary_conditions_tracer.id.end(),
ExcMessage("The boundary condition with id " +
std::to_string(triangulation_boundary_id) +
" is not present in the tracer boundary conditions."));

return (lethe_boundary_id -
this->simulation_parameters.boundary_conditions_tracer.id.begin());
}



MultiphysicsInterface<dim> *multiphysics;

TimerOutput computing_timer;
Expand Down
70 changes: 48 additions & 22 deletions source/solvers/tracer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -535,62 +535,89 @@ Tracer<dim>::assemble_system_rhs_dg()
StabilizedMethodsCopyData &copy_data)

{
// Identify which boundary condition corresponds to the boundary id. If
// this boundary condition is not identified, then exit the simulation
// instead of assuming an outlet.
const auto &triangulation_boundary_id = cell->face(face_no)->boundary_id();
const unsigned int boundary_index =
get_lethe_boundary_index(triangulation_boundary_id);
scratch_data.fe_interface_values_tracer.reinit(cell, face_no);

const FEFaceValuesBase<dim> &fe_face =
scratch_data.fe_interface_values_tracer.get_fe_face_values(0);

const auto &q_points = fe_face.get_quadrature_points();


const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
const std::vector<double> &JxW = fe_face.get_JxW_values();
const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();

std::vector<double> g(q_points.size(), 1.);

std::vector<double> values_here(q_points.size());
fe_face.get_function_values(evaluation_point, values_here);


// Gather velocity information at the face to properly advect
// First gather the dof handler for the fluid dynamics

FEFaceValues<dim> &fe_face_values_fd = scratch_data.fe_face_values_fd;

FEFaceValues<dim> &fe_face_values_fd = scratch_data.fe_face_values_fd;
const DoFHandler<dim> *dof_handler_fluid =
multiphysics->get_dof_handler(PhysicsID::fluid_dynamics);
// Identify the cell that corresponds to the fluid dynamics
typename DoFHandler<dim>::active_cell_iterator velocity_cell(
&(*triangulation), cell->level(), cell->index(), dof_handler_fluid);

fe_face_values_fd.reinit(velocity_cell, face_no);

std::vector<Tensor<1, dim>> face_velocity_values(q_points.size());
fe_face_values_fd[scratch_data.velocities].get_function_values(
*multiphysics->get_solution(PhysicsID::fluid_dynamics),
face_velocity_values);

for (unsigned int point = 0; point < q_points.size(); ++point)
{
const double velocity_dot_n =
face_velocity_values[point] * normals[point];
// If the boundary condition is an outlet, assumes that advection comes out
// since there is no inflow

if (velocity_dot_n < 0)
{
for (unsigned int i = 0; i < n_facet_dofs; ++i)
copy_data.local_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
* g[point] // g
* velocity_dot_n // \beta . n
* JxW[point]; // dx
}
if (velocity_dot_n > 0)
if (simulation_parameters.boundary_conditions_tracer.type[boundary_index] ==
BoundaryConditions::BoundaryType::outlet)
{
for (unsigned int point = 0; point < q_points.size(); ++point)
{
const double velocity_dot_n =
face_velocity_values[point] * normals[point];

for (unsigned int i = 0; i < n_facet_dofs; ++i)
copy_data.local_rhs(i) -= fe_face.shape_value(i, point) // \phi_i
* values_here[point] *
velocity_dot_n // \beta . n
* JxW[point]; // dx
}
}
// Else the boundary condition is a dirichlet. Evaluate the function and
// process it accordingly
else if (simulation_parameters.boundary_conditions_tracer
.type[boundary_index] ==
BoundaryConditions::BoundaryType::tracer_dirichlet)
{
std::vector<double> function_value(q_points.size());
simulation_parameters.boundary_conditions_tracer.tracer[boundary_index]
->value_list(q_points, function_value);


for (unsigned int point = 0; point < q_points.size(); ++point)
{
const double velocity_dot_n =
face_velocity_values[point] * normals[point];

for (unsigned int i = 0; i < n_facet_dofs; ++i)
copy_data.local_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
* function_value[point] // g
* velocity_dot_n // \beta . n
* JxW[point]; // dx
}
}
// process it accordingly
else
{
AssertThrow(false,
ExcMessage(
"No valid boundary conditions types were identified"));
}
};

const auto face_worker =
Expand Down Expand Up @@ -627,7 +654,6 @@ Tracer<dim>::assemble_system_rhs_dg()
values_there);



// Gather velocity information at the face to properly advect
// First gather the dof handler for the fluid dynamics
FEFaceValues<dim> &fe_face_values_fd = scratch_data.fe_face_values_fd;
Expand Down
47 changes: 0 additions & 47 deletions source/solvers/tracer_assemblers.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -262,31 +262,6 @@ TracerAssemblerDGCore<dim>::assemble_matrix(
}
}
} // end loop on quadrature points

// Loop over the faces to assemble the flux term
// for (unsigned int f = 0; f < scratch_data.n_faces; ++f)
// {
// for (unsigned int q = 0; q < scratch_data.n_faces_q_points; ++q)
// {
// const double flux = scratch_data.velocity_face_value[f][q] *
// scratch_data.face_normal[f][q];
// // If the flux is negative the the velocity is outbound since
// // deal.II normales point outwards In that case assemble the flux
// // if (flux > 0)
// {
// for (unsigned int i = 0; i < n_dofs; ++i)
// {
// const auto phi_T_i = scratch_data.phi_face[f][q][i];
// for (unsigned int j = 0; j < n_dofs; ++j)
// {
// const auto phi_T_j = scratch_data.phi_face[f][q][j];
// local_matrix(i, j) +=
// flux * phi_T_i * phi_T_j * scratch_data.face_JxW[f][q];
// }
// }
// }
// }
// }
}


Expand Down Expand Up @@ -331,28 +306,6 @@ TracerAssemblerDGCore<dim>::assemble_rhs(TracerScratchData<dim> &scratch_data,
JxW;
}
} // end loop on quadrature points

// Loop over the faces to assemble the flux term
// for (unsigned int f = 0; f < scratch_data.n_faces; ++f)
// {
// for (unsigned int q = 0; q < scratch_data.n_faces_q_points; ++q)
// {
// const double flux = scratch_data.velocity_face_value[f][q] *
// scratch_data.face_normal[f][q];
// // If the flux is negative the the velocity is outbound since
// // deal.II normales point outwards In that case assemble the flux
// // if (flux > 0)
// {
// for (unsigned int i = 0; i < n_dofs; ++i)
// {
// const auto phi_T_i = scratch_data.phi_face[f][q][i];
// const auto tracer_value =
// scratch_data.tracer_face_value[f][q]; local_rhs(i) -=
// flux * phi_T_i * tracer_value * scratch_data.face_JxW[f][q];
// }
// }
// }
// }
}
template class TracerAssemblerDGCore<2>;
template class TracerAssemblerDGCore<3>;
Expand Down

0 comments on commit 5e15ceb

Please sign in to comment.