Merge pull request #5481 from gassmoeller/fix_dirichlet_bcs_at_noflow…

…_boundaries

Treat noflow boundaries as inflow

doc/modules/changes/20231106_gassmoeller

@@ -0,0 +1,12 @@
+Changed: ASPECT now considers boundaries with no normal flow as
+boundaries with inflow for the purposes of the parameter
+'Allow fixed temperature on outflow boundaries' and the
+corresponding parameter for composition. This was the default
+behavior up to ASPECT 2.4.0. This behavior was changed in ASPECT 2.5.0,
+in which boundaries with no normal flow are treated like outflow
+boundaries. The new behavior caused unintended side effects, therefore
+it is reverted back to the original behavior.
+The reason for the initial change was a bugfix for boundary conditions
+in the first timestep that is now implemented differently.
+<br>
+(Rene Gassmoeller, 2023/11/06)

source/simulator/helper_functions.cc

@@ -2114,18 +2114,26 @@ namespace aspect
 
     std::vector<Tensor<1,dim>> face_current_velocity_values (fe_face_values.n_quadrature_points);
 
-    // Do not replace the id on boundaries with tangential velocity
-    const std::set<types::boundary_id> &tangential_velocity_boundaries =
+    const auto &tangential_velocity_boundaries =
       boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
 
+    const auto &zero_velocity_boundaries =
+      boundary_velocity_manager.get_zero_boundary_velocity_indicators();
+
+    const auto &prescribed_velocity_boundaries =
+      boundary_velocity_manager.get_active_boundary_velocity_conditions();
+
     // Loop over all of the boundary faces, ...
     for (const auto &cell : dof_handler.active_cell_iterators())
       if (!cell->is_artificial())
         for (const unsigned int face_number : cell->face_indices())
           {
             const typename DoFHandler<dim>::face_iterator face = cell->face(face_number);
+
+            // If the face is at a boundary where we may want to replace its id
             if (face->at_boundary() &&
-                tangential_velocity_boundaries.find(face->boundary_id()) == tangential_velocity_boundaries.end())
+                tangential_velocity_boundaries.find(face->boundary_id()) == tangential_velocity_boundaries.end() &&
+                zero_velocity_boundaries.find(face->boundary_id()) == zero_velocity_boundaries.end())
               {
                 Assert(face->boundary_id() <= offset,
                        ExcMessage("If you do not 'Allow fixed temperature/composition on outflow boundaries', "
@@ -2138,15 +2146,22 @@ namespace aspect
                 // ... check if the face is an outflow boundary by integrating the normal velocities
                 // (flux through the boundary) as: int u*n ds = Sum_q u(x_q)*n(x_q) JxW(x_q)...
                 double integrated_flow = 0;
+
                 for (unsigned int q=0; q<fe_face_values.n_quadrature_points; ++q)
                   {
-                    integrated_flow += (face_current_velocity_values[q] * fe_face_values.normal_vector(q)) *
+                    Tensor<1,dim> boundary_velocity;
+                    if (prescribed_velocity_boundaries.find(face->boundary_id()) == prescribed_velocity_boundaries.end())
+                      boundary_velocity = face_current_velocity_values[q];
+                    else
+                      boundary_velocity = boundary_velocity_manager.boundary_velocity(face->boundary_id(),
+                                                                                      fe_face_values.quadrature_point(q));
+
+                    integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
                                        fe_face_values.JxW(q);
                   }
 
                 // ... and change the boundary id of any outflow boundary faces.
-                // If there is no flow, we do not want to apply dirichlet boundary conditions either.
-                if (integrated_flow >= 0)
+                if (integrated_flow > 0)
                   face->set_boundary_id(face->boundary_id() + offset);
               }
           }

tests/airy_isostasy_initial_topo/screen-output

@@ -47,11 +47,11 @@ Number of mesh deformation degrees of freedom: 4,162
    Skipping temperature solve because RHS is zero.
    Solving C_1 system ... 9 iterations.
    Rebuilding Stokes preconditioner...
-   Solving Stokes system... 200+9 iterations.
+   Solving Stokes system... 200+10 iterations.
 
    Postprocessing:
      RMS, max velocity:    1.62e+06 m/year, 2.11e+06 m/year
-     Pressure min/avg/max: -1.118 Pa, 5.004 Pa, 10 Pa
+     Pressure min/avg/max: -1.119 Pa, 5.004 Pa, 10.01 Pa
      Topography min/max:   -4.848e-05 m, 0.04609 m
 
 *** Timestep 2:  t=3.5e-09 years, dt=1.80043e-09 years
@@ -63,7 +63,7 @@ Number of mesh deformation degrees of freedom: 4,162
 
    Postprocessing:
      RMS, max velocity:    1.53e+06 m/year, 1.99e+06 m/year
-     Pressure min/avg/max: -0.8444 Pa, 5.003 Pa, 10 Pa
+     Pressure min/avg/max: -0.8445 Pa, 5.003 Pa, 10.02 Pa
      Topography min/max:   -9.971e-05 m, 0.04229 m
 
 Termination requested by criterion: end time

tests/airy_isostasy_initial_topo/statistics

@@ -18,5 +18,5 @@
 # 18: Minimum topography (m)
 # 19: Maximum topography (m)
 0 0.000000000000e+00 0.000000000000e+00 1948 18351 8135 8135 0  0 216 798 880 1.76281656e+06 2.30007330e+06 -1.13211631e+00 5.01241125e+00 1.00000000e+01  0.00000000e+00 5.00000000e-02 
-1 1.699566636086e-09 1.699566636086e-09 1948 18351 8135 8135 0  9 209 520 860 1.62176894e+06 2.11334115e+06 -1.11840083e+00 5.00434491e+00 1.00000025e+01 -4.84780055e-05 4.60937648e-02 
-2 3.500000000000e-09 1.800433363914e-09 1948 18351 8135 8135 0 10 204 376 840 1.52731894e+06 1.98772184e+06 -8.44385424e-01 5.00258872e+00 1.00000042e+01 -9.97076582e-05 4.22912453e-02 
+1 1.699566636086e-09 1.699566636086e-09 1948 18351 8135 8135 0  9 210 580 864 1.62179316e+06 2.11339065e+06 -1.11850063e+00 5.00437021e+00 1.00114542e+01 -4.84780055e-05 4.60937648e-02 
+2 3.500000000000e-09 1.800433363914e-09 1948 18351 8135 8135 0 10 204 360 839 1.52735370e+06 1.98779365e+06 -8.44506948e-01 5.00262488e+00 1.00198524e+01 -9.97061638e-05 4.22911565e-02 

tests/artificial_viscosity_dirichlet_boundaries/screen-output

@@ -16,7 +16,7 @@ Number of degrees of freedom: 1,237 (578+81+289+289)
 
 *** Timestep 1:  t=0.00625 years, dt=0.00625 years
    Solving temperature system... 10 iterations.
-   Solving porosity system ... 11 iterations.
+   Solving porosity system ... 10 iterations.
    Rebuilding Stokes preconditioner...
    Solving Stokes system... 0+0 iterations.