Skip to content

Commit

Permalink
WIP tracer diffusivity addition
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 8349b36 commit 9fd636a
Show file tree
Hide file tree
Showing 2 changed files with 88 additions and 31 deletions.
6 changes: 3 additions & 3 deletions include/solvers/copy_data.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ class CopyData
CopyData(const unsigned int n_dofs)
: local_matrix(n_dofs, n_dofs)
, local_rhs(n_dofs)
, local_dof_indices(n_dofs) {};
, local_dof_indices(n_dofs){};

/**
* @brief Resets the cell_matrix and the cell_rhs to zero
Expand Down Expand Up @@ -96,7 +96,7 @@ class StabilizedMethodsCopyData
, local_rhs(n_dofs)
, local_dof_indices(n_dofs)
, strong_residual(n_q_points)
, strong_jacobian(n_q_points, Vector<double>(n_dofs)) {};
, strong_jacobian(n_q_points, Vector<double>(n_dofs)){};


/**
Expand Down Expand Up @@ -170,7 +170,7 @@ class StabilizedMethodsTensorCopyData
, local_rhs(n_dofs)
, local_dof_indices(n_dofs)
, strong_residual(n_q_points)
, strong_jacobian(n_q_points, std::vector<Tensor<1, dim>>(n_dofs)) {};
, strong_jacobian(n_q_points, std::vector<Tensor<1, dim>>(n_dofs)){};

/**
* @brief Resets the cell_matrix, cell_rhs, strong_residual
Expand Down
113 changes: 85 additions & 28 deletions source/solvers/tracer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@ Tracer<dim>::assemble_system_matrix_dg()
const unsigned int &face_no,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
const double beta = 10;
const double beta = 0;
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);
Expand All @@ -245,6 +245,18 @@ Tracer<dim>::assemble_system_matrix_dg()
const std::vector<double> &JxW = fe_face.get_JxW_values();
const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();



// Calculate diffusivity at the faces
auto &properties_manager =
this->simulation_parameters.physical_properties_manager;
const auto diffusivity_model =
properties_manager.get_tracer_diffusivity();
std::map<field, std::vector<double>> fields;
std::vector<double> tracer_diffusivity(q_points.size());
diffusivity_model->vector_value(fields, tracer_diffusivity);


// 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 All @@ -267,25 +279,28 @@ Tracer<dim>::assemble_system_matrix_dg()
const double velocity_dot_n =
face_velocity_values[point] * normals[point];

if (velocity_dot_n > 0)
{
for (unsigned int i = 0; i < n_facet_dofs; ++i)
for (unsigned int j = 0; j < n_facet_dofs; ++j)

for (unsigned int i = 0; i < n_facet_dofs; ++i)
for (unsigned int j = 0; j < n_facet_dofs; ++j)
{
if (velocity_dot_n > 0)
{
copy_data.local_matrix(i, j) +=
fe_face.shape_value(i, point) // \phi_i
* fe_face.shape_value(j, point) // \phi_j
* velocity_dot_n // \beta . n
* JxW[point]; // dx
}


copy_data.local_matrix(i, j) +=
fe_face.shape_value(i, point) // \phi_i
* fe_face.shape_value(j, point) // \phi_j
* beta // \beta . n
* JxW[point]; // dx
}
}
copy_data.local_matrix(i, j) -=
tracer_diffusivity[point] *
(fe_face.shape_value(i, point) *
fe_face.shape_grad(j, point) * normals[point] +
fe_face.shape_value(j, point) *
fe_face.shape_grad(i, point) * normals[point]) *
JxW[point];
}
}
};

Expand All @@ -298,7 +313,7 @@ Tracer<dim>::assemble_system_matrix_dg()
const unsigned int &nsf,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
const double beta = 10;
const double beta = 0;
FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values_tracer;
fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
const auto &q_points = fe_iv.get_quadrature_points();
Expand All @@ -314,6 +329,14 @@ Tracer<dim>::assemble_system_matrix_dg()
const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();


// Calculate diffusivity at the faces
auto &properties_manager =
this->simulation_parameters.physical_properties_manager;
const auto diffusivity_model =
properties_manager.get_tracer_diffusivity();
std::map<field, std::vector<double>> fields;
std::vector<double> tracer_diffusivity(q_points.size());
diffusivity_model->vector_value(fields, tracer_diffusivity);

// Gather velocity information at the face to properly advect
// First gather the dof handler for the fluid dynamics
Expand Down Expand Up @@ -349,8 +372,12 @@ Tracer<dim>::assemble_system_matrix_dg()
// symmetric interior penalty method. See Larson Chap. 14. P.362

copy_data_face.cell_matrix(i, j) +=
(-fe_iv.average_of_shape_gradients(j, q) * normals[q] *
fe_iv.jump_in_shape_values(i, q) +
(-tracer_diffusivity[q] *
fe_iv.average_of_shape_gradients(j, q) * normals[q] *
fe_iv.jump_in_shape_values(i, q) -
tracer_diffusivity[q] *
fe_iv.average_of_shape_gradients(i, q) * normals[q] *
fe_iv.jump_in_shape_values(j, q) +
beta * fe_iv.jump_in_shape_values(j, q) *
fe_iv.jump_in_shape_values(i, q)) *
JxW[q];
Expand Down Expand Up @@ -559,6 +586,8 @@ Tracer<dim>::assemble_system_rhs_dg()
StabilizedMethodsCopyData &copy_data)

{
const double beta = 0;

// 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.
Expand All @@ -580,6 +609,14 @@ Tracer<dim>::assemble_system_rhs_dg()
std::vector<double> values_here(q_points.size());
fe_face.get_function_values(evaluation_point, values_here);

// Calculate diffusivity at the faces
auto &properties_manager =
this->simulation_parameters.physical_properties_manager;
const auto diffusivity_model = properties_manager.get_tracer_diffusivity();
std::map<field, std::vector<double>> fields;
std::vector<double> tracer_diffusivity(q_points.size());
diffusivity_model->vector_value(fields, tracer_diffusivity);

// 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 Expand Up @@ -630,11 +667,17 @@ Tracer<dim>::assemble_system_rhs_dg()

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
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

// copy_data.local_rhs(i) -=
// fe_face.shape_value(i, point) // \phi_i
// * function_value[point] // \phi_j
// * beta // \beta . n
// * JxW[point]; // dx
}
}
}
Expand All @@ -657,11 +700,12 @@ Tracer<dim>::assemble_system_rhs_dg()
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
// TODO refactor and put inside a parameter formally
const double beta = 10.;
const double beta = 0.;
FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values_tracer;
fe_iv.reinit(cell, f, sf, ncell, nf, nsf);



const auto &q_points = fe_iv.get_quadrature_points();

// Allocate memory for the faces calculation
Expand All @@ -685,7 +729,18 @@ Tracer<dim>::assemble_system_rhs_dg()
std::vector<Tensor<1, dim>> tracer_average_gradient(q_points.size());
std::vector<double> tracer_value_jump(q_points.size());

fe_iv.get_average_of_function_values(evaluation_point, tracer_value_jump);
// Calculate diffusivity at the faces
auto &properties_manager =
this->simulation_parameters.physical_properties_manager;
const auto diffusivity_model =
properties_manager.get_tracer_diffusivity();
std::map<field, std::vector<double>> fields;
std::vector<double> tracer_diffusivity(q_points.size());
diffusivity_model->vector_value(fields, tracer_diffusivity);


fe_iv.get_jump_in_function_values(evaluation_point, tracer_value_jump);

fe_iv.get_average_of_function_gradients(evaluation_point,
tracer_average_gradient);

Expand Down Expand Up @@ -734,8 +789,10 @@ Tracer<dim>::assemble_system_rhs_dg()
// penalty method. See Larson Chap. 14. P.362

copy_data_face.cell_rhs(i) -=
(-tracer_average_gradient[q] * normals[q] *
fe_iv.jump_in_shape_values(i, q) +
(-tracer_diffusivity[q] * tracer_average_gradient[q] *
normals[q] * fe_iv.jump_in_shape_values(i, q) -
tracer_diffusivity[q] * tracer_value_jump[q] * normals[q] *
fe_iv.average_of_shape_gradients(i, q) +
beta * tracer_value_jump[q] *
fe_iv.jump_in_shape_values(i, q)) *
JxW[q];
Expand Down Expand Up @@ -1148,10 +1205,10 @@ Tracer<dim>::postprocess_tracer_flow_rate(const VectorType &current_solution_fd)
normal_vector_tracer) *
fe_face_values_tracer.JxW(q);
} // end loop on quadrature points
} // end face is a boundary face
} // end loop on faces
} // end condition cell at boundary
} // end loop on cells
} // end face is a boundary face
} // end loop on faces
} // end condition cell at boundary
} // end loop on cells


// Sum across all cores
Expand Down

0 comments on commit 9fd636a

Please sign in to comment.