Skip to content

Commit

Permalink
WIP Yo
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 5e15ceb commit 8349b36
Showing 1 changed file with 77 additions and 30 deletions.
107 changes: 77 additions & 30 deletions source/solvers/tracer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,7 @@ Tracer<dim>::assemble_system_matrix_dg()
const unsigned int &face_no,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
const double beta = 10;
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 Down Expand Up @@ -270,11 +271,20 @@ Tracer<dim>::assemble_system_matrix_dg()
{
for (unsigned int i = 0; i < n_facet_dofs; ++i)
for (unsigned int j = 0; j < n_facet_dofs; ++j)
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
* 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
}
}
}
};
Expand All @@ -288,6 +298,7 @@ Tracer<dim>::assemble_system_matrix_dg()
const unsigned int &nsf,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
const double beta = 10;
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 @@ -302,6 +313,8 @@ Tracer<dim>::assemble_system_matrix_dg()
const std::vector<double> &JxW = fe_iv.get_JxW_values();
const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();



// 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 @@ -318,19 +331,30 @@ Tracer<dim>::assemble_system_matrix_dg()
*multiphysics->get_solution(PhysicsID::fluid_dynamics),
face_velocity_values);

for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
for (unsigned int q = 0; q < q_points.size(); ++q)
{
const double velocity_dot_n =
face_velocity_values[qpoint] * normals[qpoint];
const double velocity_dot_n = face_velocity_values[q] * normals[q];
for (unsigned int i = 0; i < n_dofs; ++i)
for (unsigned int j = 0; j < n_dofs; ++j)
copy_data_face.cell_matrix(i, j) +=
fe_iv.jump_in_shape_values(i, qpoint) // [\phi_i]
* fe_iv.shape_value((velocity_dot_n > 0),
j,
qpoint) // phi_j^{upwind}
* velocity_dot_n // (\beta .n)
* JxW[qpoint]; // dx
{
copy_data_face.cell_matrix(i, j) +=
fe_iv.jump_in_shape_values(i, q) // [\phi_i]
* fe_iv.shape_value((velocity_dot_n > 0),
j,
q) // phi_j^{upwind}
* velocity_dot_n // (\beta .n)
* JxW[q]; // dx

// Assemble the diffusion term using nitsche
// 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) +
beta * fe_iv.jump_in_shape_values(j, q) *
fe_iv.jump_in_shape_values(i, q)) *
JxW[q];
}
}
};

Expand Down Expand Up @@ -605,10 +629,13 @@ Tracer<dim>::assemble_system_rhs_dg()
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] // g
* velocity_dot_n // \beta . n
* JxW[point]; // dx
}
}
}
// process it accordingly
Expand All @@ -629,6 +656,8 @@ Tracer<dim>::assemble_system_rhs_dg()
const unsigned int &nsf,
TracerScratchData<dim> &scratch_data,
StabilizedMethodsCopyData &copy_data) {
// TODO refactor and put inside a parameter formally
const double beta = 10.;
FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values_tracer;
fe_iv.reinit(cell, f, sf, ncell, nf, nsf);

Expand All @@ -653,6 +682,14 @@ Tracer<dim>::assemble_system_rhs_dg()
fe_iv.get_fe_face_values(1).get_function_values(evaluation_point,
values_there);

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);
fe_iv.get_average_of_function_gradients(evaluation_point,
tracer_average_gradient);



// Gather velocity information at the face to properly advect
// First gather the dof handler for the fluid dynamics
Expand All @@ -670,28 +707,38 @@ Tracer<dim>::assemble_system_rhs_dg()
*multiphysics->get_solution(PhysicsID::fluid_dynamics),
face_velocity_values);

for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
for (unsigned int q = 0; q < q_points.size(); ++q)
{
const double velocity_dot_n =
face_velocity_values[qpoint] * normals[qpoint];
const double velocity_dot_n = face_velocity_values[q] * normals[q];
for (unsigned int i = 0; i < n_dofs; ++i)
{
// Assemble advection terms with upwinding
if (velocity_dot_n > 0)
{
copy_data_face.cell_rhs(i) -=
fe_iv.jump_in_shape_values(i, qpoint) // [\phi_i]
* values_here[qpoint] // \phi^{upwind}
* velocity_dot_n // (\beta .n)
* JxW[qpoint]; // dx
fe_iv.jump_in_shape_values(i, q) // [\phi_i]
* values_here[q] // \phi^{upwind}
* velocity_dot_n // (\beta .n)
* JxW[q]; // dx
}
else
{
copy_data_face.cell_rhs(i) -=
fe_iv.jump_in_shape_values(i, qpoint) // [\phi_i]
* values_there[qpoint] // \phi^{upwind}
* velocity_dot_n // (\beta .n)
* JxW[qpoint]; // dx
fe_iv.jump_in_shape_values(i, q) // [\phi_i]
* values_there[q] // \phi^{upwind}
* velocity_dot_n // (\beta .n)
* JxW[q]; // dx
}

// Assemble the diffusion term using nitsche symmetric interior
// 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) +
beta * tracer_value_jump[q] *
fe_iv.jump_in_shape_values(i, q)) *
JxW[q];
}
}
};
Expand Down

0 comments on commit 8349b36

Please sign in to comment.