A Robin boundary condition for structure module #555
Conversation
| std::map<dealii::types::boundary_id, std::pair<std::array<bool, 2>, std::array<double, 3>>> | ||
| robin_k_c_p_param; |
There was a problem hiding this comment.
This design makes the code very hard to read in my opinion, because in the implementation (SpatialOperator classes), one only sees .first, .second, etc., without concrete names.
Here, the data structure should look like
std::map<dealii::types::boundary_id, ExaDG::RobinParameters> robin_bc;and inside RobinParameters variables like .stiffness, etc. with readable names. Instead of robin_k_c_p_param (where it is hard to guess the meaning), we should follow the naming robin_bc for consistency with other boundary types.
| // integrate standard (stored) traction or exterior pressure on Robin boundaries | ||
| if(boundary_type == BoundaryType::Neumann or boundary_type == BoundaryType::NeumannCached or | ||
| boundary_type == BoundaryType::RobinSpringDashpotPressure) | ||
| { | ||
| tensor F = get_F<dim, Number>(integrator.get_gradient(q)); | ||
| vector N = integrator.get_normal_vector(q); | ||
| // da/dA * n = det F F^{-T} * N := n_star | ||
| // -> da/dA = n_star.norm() | ||
| vector n_star = determinant(F) * transpose(invert(F)) * N; | ||
| // t_0 = da/dA * t | ||
| traction *= n_star.norm(); | ||
| if(operator_type == OperatorType::inhomogeneous or operator_type == OperatorType::full) | ||
| { |
There was a problem hiding this comment.
From the many if-statements introduced in this function, I am not sure whether we have a good design right now for the new Robin-functionality in the present PR. Note that the free function calculate_neumann_value() already contains if-else-statements. I am also not convinced right now that the way OperatorBase is extended is a good design.
| traction *= n_star.norm(); | ||
| if(operator_type == OperatorType::inhomogeneous or operator_type == OperatorType::full) | ||
| { | ||
| traction -= calculate_neumann_value<dim, Number>( |
There was a problem hiding this comment.
I don't like that a minus sign is introduced here, and that the sign is also changed below in function submit_value().
(As a sidenote, the present change turns an assignment into an addition/subtraction, so why didn't you actually write traction = - xyz; to introduce a minus sign compared to the previous version? However, we should keep the previous version in my opinion. Robin is a new boundary type, that is optional.)
|
|
||
| if(this->operator_data.pull_back_traction) | ||
| // integrate standard (stored) traction or exterior pressure on Robin boundaries |
There was a problem hiding this comment.
This comment reads like Robin is "the only" boundary type. I think it holds for the present PR as a whole: We should aim at a modular code where Robin is one specific feature and where it would be very easy in the future to remove this feature again (with "easy" I mean minimal code changes). Even if we will never remove it again, this way of thinking might help for the design.
| if(param.problem_type == ProblemType::Unsteady) | ||
| { | ||
| // The boundary mass operator is only required, if we have damping terms from the Robin | ||
| // condition, i.e., in the unsteady case. | ||
| matrix_free_data.append_mapping_flags( | ||
| BoundaryMassOperator<dim, Number, dim /* n_components */>::get_mapping_flags()); | ||
| } |
There was a problem hiding this comment.
Similar to this comment https://github.com/exadg/exadg/pull/555/files#r1457203165: I think we do not always have Robin type BCs, so I would expect an if-statement here like if(robin){...}?
| // check boundary ID in robin_k_c_p_param to add boundary mass integrals from Robin boundaries | ||
| // on BoundaryType::NeumannCached or BoundaryType::RobinSpringDashpotPressure | ||
| if(boundary_type == BoundaryType::NeumannCached or | ||
| boundary_type == BoundaryType::RobinSpringDashpotPressure) | ||
| { | ||
| if(operator_type == OperatorType::homogeneous or operator_type == OperatorType::full) | ||
| { | ||
| auto const it = this->operator_data.bc->robin_k_c_p_param.find(boundary_id); | ||
|
|
||
| if(it != this->operator_data.bc->robin_k_c_p_param.end()) | ||
| { | ||
| bool const normal_projection_displacement = it->second.first[0]; | ||
| double const coefficient_displacement = it->second.second[0]; | ||
|
|
||
| if(normal_projection_displacement) | ||
| { | ||
| vector const N = integrator.get_normal_vector(q); | ||
| traction += N * (coefficient_displacement * (N * integrator.get_value(q))); | ||
| } | ||
| else | ||
| { | ||
| traction += coefficient_displacement * integrator.get_value(q); | ||
| } | ||
|
|
||
| if(this->operator_data.unsteady) | ||
| { | ||
| bool const normal_projection_velocity = it->second.first[1]; | ||
| double const coefficient_velocity = it->second.second[1]; | ||
|
|
||
| if(normal_projection_velocity) | ||
| { | ||
| vector const N = integrator.get_normal_vector(q); | ||
| traction += N * (coefficient_velocity * this->scaling_factor_mass_boundary * | ||
| (N * integrator.get_value(q))); | ||
| } | ||
| else | ||
| { | ||
| traction += | ||
| coefficient_velocity * this->scaling_factor_mass_boundary * integrator.get_value(q); | ||
| } | ||
| } | ||
| } | ||
| } |
There was a problem hiding this comment.
This code looks like a duplication of the one above in linear_operator.cpp. Is there any difference between these two version?
In the previous code version, one could say we had a function call calculate_neumann_value() in both cases (linear operator and nonlinear operator), but here we have a sequence of instructions and calculations that are duplicated, which should definitely be avoided.
| dealii::AlignedVector<dealii::VectorizedArray<Number>> local_diag(dofs_per_cell); | ||
|
|
||
| for(auto face = range.first; face < range.second; ++face) | ||
| if(is_dg) |
There was a problem hiding this comment.
I think this is a mixture of different styles. In other places, we use functions face_loop_empty() so that we see on the top level (matrix_free->loop()/cell_loop()) which type of loops we have. I think we need a consistent way, otherwise we can here and there introduce if-statements "for safety". So I think we should not write if(is_dg) here.
| @@ -823,10 +861,10 @@ OperatorBase<dim, Number, n_components>::internal_calculate_system_matrix( | |||
| SparseMatrix & system_matrix) const | |||
| { | |||
| // assemble matrix locally on each process | |||
| if(evaluate_face_integrals() and is_dg) | |||
| if(evaluate_face_integrals()) | |||
There was a problem hiding this comment.
I find this logic irritating or not consistent with other parts of OperatorBase where we pass empty functions to matrix_free->loop(). I think we should distinguish between dg and cg on the top level, and inside each variant introduce another if depending on whether we have face loops or not. If we enter a face_loop, we will then also do the face integrals. See also comment https://github.com/exadg/exadg/pull/555/files#r1457261288
|
I think it makes sense to split this PR into parts as you suggest above. The starting point are probably changes made to From the code review I did so far, I see mainly three PRs (1. OperatorBase, 2. introduction of BoundaryMassOperator, 3. all the rest; but we should do them one after the other and then the picture might also become clearer). |
This PR introduces a general BC of the form
traction = P*N = - k * displacement - c * velocity - exterior_pressure.For this we need several adaptations:
boundary_mass_operatorto shift the discretized remainder of the velocity term to the rhs (similar to the weak damping or general time integration terms, but now with a mass term on the boundary only). This will also come in handy, if we want to project values on the boundary.operator_basesuch that CG discretizations involve some boundary integrals ifevaluate_boundary_integrals() == truein the fundamentalapply,apply_add,evaluate,evaluate_add... operations and the remaining MF loops used in MG setup, preconditioner update etc. (building the system matrix or its inverse diagonal)BoundaryMassOperatorinstructure/operator.cppand the according interface.I have this now in a single PR to thoroughly check if the behavior is as desired (see the tests) and the changes in
OperatorBaseactually work as expected to get good Multigrid convergence. Also, with this you can see the intended use (and extensibility towards Robin coupling conditions in FSI). If need be, I can split into several parts:OperatorBaseBoundaryMassOperatorBoundaryMassOperatorand temporal scaling)Whoever is reviewing this (probably @nfehn 馃), let's have an in-person discussion first!