Skip to content

Commit

Permalink
WIP smthing weird
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 9fd636a commit 0aab39a
Showing 1 changed file with 134 additions and 119 deletions.
253 changes: 134 additions & 119 deletions source/solvers/tracer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,10 @@ Tracer<dim>::assemble_system_matrix_dg()
const unsigned int &face_no,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
const double beta = 0;
double beta = 1;
beta *= 1 / compute_cell_diameter<dim>(cell->measure(), 1);


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,8 +248,6 @@ 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;
Expand Down Expand Up @@ -286,20 +287,21 @@ Tracer<dim>::assemble_system_matrix_dg()
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
fe_face.shape_value(i, point) *
fe_face.shape_value(j, point) * velocity_dot_n *
JxW[point];
}


copy_data.local_matrix(i, j) -=
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];
(-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] +
beta * fe_face.shape_value(i, point) *
fe_face.shape_value(j, point) * JxW[point];
}
}
};
Expand All @@ -313,7 +315,8 @@ Tracer<dim>::assemble_system_matrix_dg()
const unsigned int &nsf,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
const double beta = 0;
double beta = 0;
beta *= 1 / compute_cell_diameter<dim>(cell->measure(), 1);
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 Down Expand Up @@ -583,112 +586,124 @@ Tracer<dim>::assemble_system_rhs_dg()
[&](const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int &face_no,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data)
StabilizedMethodsCopyData &copy_data) {
double beta = 1;

{
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.
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> 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;
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);

// If the boundary condition is an outlet, assumes that advection comes out
// since there is no inflow

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);
beta *= 1 / compute_cell_diameter<dim>(cell->measure(), 1);


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)
{
// 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> values_here(q_points.size());
std::vector<Tensor<1, dim>> gradients_here(q_points.size());

fe_face.get_function_values(evaluation_point, values_here);
fe_face.get_function_gradients(evaluation_point, gradients_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;
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);

// If the boundary condition is an outlet, assumes that advection comes
// out since there is no inflow

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
* 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
}
}
}
// process it accordingly
else
{
AssertThrow(false,
ExcMessage(
"No valid boundary conditions types were identified"));
}
};
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) *
function_value[point] *
velocity_dot_n * JxW[point];

copy_data.local_rhs(i) -=
tracer_diffusivity[point] *
(-fe_face.shape_value(i, point) * gradients_here[point] *
normals[point] -
(values_here[point] - function_value[point]) *
fe_face.shape_grad(i, point) * normals[point]) *
JxW[point] +
beta * fe_face.shape_value(i, point) *
(values_here[point] - function_value[point]) * JxW[point];
}
}
}
// process it accordingly
else
{
AssertThrow(false,
ExcMessage(
"No valid boundary conditions types were identified"));
}
};

const auto face_worker =
[&](const typename DoFHandler<dim>::active_cell_iterator &cell,
Expand Down Expand Up @@ -1205,10 +1220,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 0aab39a

Please sign in to comment.