diff --git a/source/solvers/tracer.cc b/source/solvers/tracer.cc index 652b99222c..b7e1c4031a 100644 --- a/source/solvers/tracer.cc +++ b/source/solvers/tracer.cc @@ -234,7 +234,10 @@ Tracer::assemble_system_matrix_dg() const unsigned int &face_no, TracerScratchData &scratch_data, StabilizedMethodsCopyData ©_data) { - const double beta = 0; + double beta = 1; + beta *= 1 / compute_cell_diameter(cell->measure(), 1); + + scratch_data.fe_interface_values_tracer.reinit(cell, face_no); const FEFaceValuesBase &fe_face = scratch_data.fe_interface_values_tracer.get_fe_face_values(0); @@ -245,8 +248,6 @@ Tracer::assemble_system_matrix_dg() const std::vector &JxW = fe_face.get_JxW_values(); const std::vector> &normals = fe_face.get_normal_vectors(); - - // Calculate diffusivity at the faces auto &properties_manager = this->simulation_parameters.physical_properties_manager; @@ -286,20 +287,21 @@ Tracer::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]; } } }; @@ -313,7 +315,8 @@ Tracer::assemble_system_matrix_dg() const unsigned int &nsf, TracerScratchData &scratch_data, StabilizedMethodsCopyData ©_data) { - const double beta = 0; + double beta = 0; + beta *= 1 / compute_cell_diameter(cell->measure(), 1); FEInterfaceValues &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(); @@ -583,112 +586,124 @@ Tracer::assemble_system_rhs_dg() [&](const typename DoFHandler::active_cell_iterator &cell, const unsigned int &face_no, TracerScratchData &scratch_data, - StabilizedMethodsCopyData ©_data) + StabilizedMethodsCopyData ©_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 &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 &JxW = fe_face.get_JxW_values(); - const std::vector> &normals = fe_face.get_normal_vectors(); - - std::vector 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> fields; - std::vector 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 &fe_face_values_fd = scratch_data.fe_face_values_fd; - const DoFHandler *dof_handler_fluid = - multiphysics->get_dof_handler(PhysicsID::fluid_dynamics); - // Identify the cell that corresponds to the fluid dynamics - typename DoFHandler::active_cell_iterator velocity_cell( - &(*triangulation), cell->level(), cell->index(), dof_handler_fluid); - fe_face_values_fd.reinit(velocity_cell, face_no); - std::vector> 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 function_value(q_points.size()); - simulation_parameters.boundary_conditions_tracer.tracer[boundary_index] - ->value_list(q_points, function_value); + beta *= 1 / compute_cell_diameter(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 &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 &JxW = fe_face.get_JxW_values(); + const std::vector> &normals = fe_face.get_normal_vectors(); + + std::vector values_here(q_points.size()); + std::vector> 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> fields; + std::vector 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 &fe_face_values_fd = scratch_data.fe_face_values_fd; + const DoFHandler *dof_handler_fluid = + multiphysics->get_dof_handler(PhysicsID::fluid_dynamics); + // Identify the cell that corresponds to the fluid dynamics + typename DoFHandler::active_cell_iterator velocity_cell( + &(*triangulation), cell->level(), cell->index(), dof_handler_fluid); + fe_face_values_fd.reinit(velocity_cell, face_no); + std::vector> 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 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::active_cell_iterator &cell, @@ -1205,10 +1220,10 @@ Tracer::postprocess_tracer_flow_rate(const VectorType ¤t_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