Merge pull request #6067 from anne-glerum/adjust_v_for_outflow_BC

Subtract mesh velocity before computing whether a boundary experiences in- or outflow for fixed boundary conditions

benchmarks/diffusion_of_hill/1_sine_zero_flux_ci.prm

@@ -0,0 +1,162 @@
+set Additional shared libraries            = libanalytical_topography.so
+set Dimension                              = 2
+set Use years in output instead of seconds = false
+set End time                               = 0.005
+set Maximum time step                      = 0.0005
+set Output directory                       = output-1_sine_zero_flux_ci_toppresvel_subtractmeshvel/
+set Nonlinear solver scheme = single Advection, single Stokes
+set Pressure normalization                 = surface
+set Surface pressure                       = 0
+
+# 1x1 box with an initial hill topography
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+  end
+
+  subsection Initial topography model
+    set Model name = function
+
+    subsection Function
+      set Function constants = A=0.150, L=1.
+      set Function expression = \
+                                if(x<0.5,A * sin((x+0.5)*pi),0)
+    end
+  end
+end
+
+# Temperature effects are ignored
+subsection Initial temperature model
+  set Model name = function
+
+  subsection Function
+    set Function expression = 0
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = bottom, top
+  set List of model names = initial temperature
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+  set Names of fields = layer_1, layer_2, layer_3
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names = x,z
+    set Coordinate system = cartesian #depth
+    #set Function expression = if(x<.05, 1, 0); if(x>=0.05 && x<0.1, 1, 0); 0
+    set Function expression = if(z>1.07, 1, 0); 0; 0
+  end
+end
+
+subsection Boundary composition model
+  set Fixed composition boundary indicators = top
+  set Allow fixed composition on outflow boundaries = false
+  set List of model names = function
+  subsection Function
+    set Variable names = x,z,t
+    set Function expression = if(t==0,1,0); 0; if(t==0,0,1)
+  end
+end
+
+# Free slip on all boundaries
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, bottom
+  set Prescribed velocity boundary indicators = top: function
+  subsection Function
+    set Variable names = x,z,t
+    set Function expression = 0;0
+  end
+end
+
+# The mesh will be deformed according to the displacement
+# of the surface nodes due to diffusion of the topography.
+# The mesh is allowed to move vertical along the left and
+# right boundary.
+subsection Mesh deformation
+  set Mesh deformation boundary indicators = top: diffusion
+  set Additional tangential mesh velocity boundary indicators = left, right
+
+  subsection Diffusion
+    # The diffusivity
+    set Hillslope transport coefficient = 0.25
+  end
+end
+
+# Vertical gravity
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0
+  end
+end
+
+# One material with unity properties
+subsection Material model
+  set Model name = simple
+
+  subsection Simple model
+    set Reference density             = 1
+    set Reference specific heat       = 1
+    set Reference temperature         = 0
+    set Thermal conductivity          = 1
+    set Thermal expansion coefficient = 1
+    set Viscosity                     = 1
+  end
+end
+
+# We also have to specify that we want to use the Boussinesq
+# approximation (assuming the density in the temperature
+# equation to be constant, and incompressibility).
+subsection Formulation
+  set Formulation = Boussinesq approximation
+end
+
+# We here use a globally refined mesh without
+# adaptive mesh refinement.
+subsection Mesh refinement
+  set Initial global refinement                = 6
+  set Initial adaptive refinement              = 0
+  set Time steps between mesh refinement       = 0
+  set Strategy = minimum refinement function
+  subsection Minimum refinement function
+    set Variable names = x,z
+    set Coordinate system = cartesian
+    set Function expression = if(z>0.87,7,3)
+  end
+end
+
+# We output the computed topography and the analytical topography
+# value to file.
+subsection Postprocess
+  set List of postprocessors = velocity statistics, temperature statistics, heat flux statistics, visualization #, analytical topography
+
+  subsection Topography
+    set Output to file = true
+    set Analytical solution of example = 1
+    set Diffusivity = 0.25
+    set Initial sinusoidal topography amplitude = 0.075
+  end
+
+  subsection Visualization
+    set Time between graphical output = 0.0001
+    set Output mesh velocity = true
+    set Interpolate output = false
+  end
+end
+
+subsection Solver parameters
+  subsection Stokes solver parameters
+    set Linear solver tolerance                         = 1e-7
+  end
+end

doc/modules/changes/20241111_glerum

@@ -0,0 +1,8 @@
+ Fixed: The function to determine whether a boundary is an inflow
+or outflow boundary (for the purpose of deciding whether to apply
+boundary conditions for compositional fields) did not take into
+account the mesh deformation. This caused incorrectly applied
+surface boundary conditions in models with free surface and is
+fixed now.
+ <br>
+ (Anne Glerum, 2024/11/11)

source/simulator/helper_functions.cc

@@ -2306,6 +2306,7 @@ namespace aspect
                                       update_quadrature_points | update_JxW_values);
 
     std::vector<Tensor<1,dim>> face_current_velocity_values (fe_face_values.n_quadrature_points);
+    std::vector<Tensor<1,dim>> face_current_mesh_velocity_values (fe_face_values.n_quadrature_points);
 
     const auto &tangential_velocity_boundaries =
       boundary_velocity_manager.get_tangential_boundary_velocity_indicators();
@@ -2335,6 +2336,10 @@ namespace aspect
                 fe_face_values.reinit (cell, face_number);
                 fe_face_values[introspection.extractors.velocities].get_function_values(current_linearization_point,
                                                                                         face_current_velocity_values);
+                // get the mesh velocity, as we need to subtract it off of the advection systems
+                if (parameters.mesh_deformation_enabled)
+                  fe_face_values[introspection.extractors.velocities].get_function_values(mesh_deformation->mesh_velocity,
+                                                                                          face_current_mesh_velocity_values);
 
                 // ... 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)...
@@ -2349,6 +2354,9 @@ namespace aspect
                       boundary_velocity = boundary_velocity_manager.boundary_velocity(face->boundary_id(),
                                                                                       fe_face_values.quadrature_point(q));
 
+                    if (parameters.mesh_deformation_enabled)
+                      boundary_velocity -= face_current_mesh_velocity_values[q];
+
                     integrated_flow += (boundary_velocity * fe_face_values.normal_vector(q)) *
                                        fe_face_values.JxW(q);
                   }

tests/hill_diffusion_field_boundary_conditions.prm

@@ -0,0 +1,98 @@
+# This setup tests the compositional field boundary conditions
+# along a deforming boundary. To determine whether a part of the
+# boundary experiences in- or outflow, the mesh velocity is
+# subtracted from the material velocity. In this test, the
+# prescribed material velocity is zero, and only the erosion
+# or deposition (decrease or increase in height) by the diffusion
+# applied to the surface topography determines whether boundary
+# conditions should be applied (inflow/deposition),
+# or not (outflow/erosion).
+
+include $ASPECT_SOURCE_DIR/tests/hill_diffusion.prm
+
+set Dimension                              = 2
+set Nonlinear solver scheme                = single Advection, single Stokes
+
+# 1x1 box with an initial hill topography
+subsection Geometry model
+  set Model name = box
+
+  subsection Box
+    set X extent = 1
+    set Y extent = 1
+  end
+
+  subsection Initial topography model
+    set Model name = function
+
+    subsection Function
+      set Function constants = A=0.150, L=1.
+      set Function expression = \
+                                if(x<0.5,A * sin((x+0.5)*pi),0)
+    end
+  end
+end
+
+# We here use a globally refined mesh without
+# adaptive mesh refinement.
+subsection Mesh refinement
+  set Initial global refinement                = 4
+end
+
+subsection Compositional fields
+  set Number of fields = 3
+  set Names of fields = layer_1, layer_2, layer_3
+end
+
+subsection Initial composition model
+  set Model name = function
+
+  subsection Function
+    set Variable names = x,z
+    set Coordinate system = cartesian #depth
+    #set Function expression = if(x<.05, 1, 0); if(x>=0.05 && x<0.1, 1, 0); 0
+    set Function expression = if(z>1.07, 1, 0); 0; 0
+  end
+end
+
+subsection Boundary composition model
+  set Fixed composition boundary indicators = top
+  set Allow fixed composition on outflow boundaries = false
+  set List of model names = function
+  subsection Function
+    set Variable names = x,z,t
+    set Function expression = if(t==0,1,0); 0; if(t==0,0,1)
+  end
+end
+
+# Free slip on all boundaries except the top boundary.
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, bottom
+  set Prescribed velocity boundary indicators = top: function
+  subsection Function
+    set Variable names = x,z,t
+    set Function expression = 0;0
+  end
+end
+
+# Vertical gravity
+subsection Gravity model
+  set Model name = vertical
+
+  subsection Vertical
+    set Magnitude = 0
+  end
+end
+
+# We output the computed topography
+# value to file.
+subsection Postprocess
+  set List of postprocessors = velocity statistics, topography, visualization
+
+  subsection Visualization
+    set Time between graphical output = 0.0001
+    set Output mesh velocity = true
+    set Interpolate output = false
+    set Output format = gnuplot
+  end
+end

tests/hill_diffusion_field_boundary_conditions/screen-output

@@ -0,0 +1,63 @@
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------
+Number of active cells: 256 (on 5 levels)
+Number of degrees of freedom: 6,823 (2,178+289+1,089+1,089+1,089+1,089)
+
+Number of mesh deformation degrees of freedom: 578
+   Solving mesh displacement system... 0 iterations.
+   Added initial topography to grid
+
+*** Timestep 0:  t=0 seconds, dt=0 seconds
+   Solving mesh displacement system... 0 iterations.
+   Skipping temperature solve because RHS is zero.
+   Solving layer_1 system ... 0 iterations.
+   Skipping layer_2 composition solve because RHS is zero.
+   Skipping layer_3 composition solve because RHS is zero.
+   Solving Stokes system (GMG)... 0+0 iterations.
+
+   Postprocessing:
+     RMS, max velocity:        0 m/s, 0 m/s
+     Topography min/max:       0 m, 0.15 m
+     Writing graphical output: output-hill_diffusion_field_boundary_conditions/solution/solution-00000
+
+*** Timestep 1:  t=0.0005 seconds, dt=0.0005 seconds
+   Solving mesh surface diffusion
+   Solving mesh displacement system... 4 iterations.
+   Skipping temperature solve because RHS is zero.
+   Solving layer_1 system ... 11 iterations.
+   Skipping layer_2 composition solve because RHS is zero.
+   Solving layer_3 system ... 12 iterations.
+   Solving Stokes system (GMG)... 0+0 iterations.
+
+   Postprocessing:
+     RMS, max velocity:        0 m/s, 0 m/s
+     Topography min/max:       -0.0002308 m, 0.1498 m
+     Writing graphical output: output-hill_diffusion_field_boundary_conditions/solution/solution-00001
+
+*** Timestep 2:  t=0.001 seconds, dt=0.0005 seconds
+   Solving mesh surface diffusion
+   Solving mesh displacement system... 4 iterations.
+   Skipping temperature solve because RHS is zero.
+   Solving layer_1 system ... 12 iterations.
+   Skipping layer_2 composition solve because RHS is zero.
+   Solving layer_3 system ... 12 iterations.
+   Solving Stokes system (GMG)... 0+0 iterations.
+
+   Postprocessing:
+     RMS, max velocity:        0 m/s, 0 m/s
+     Topography min/max:       -0.0003541 m, 0.1496 m
+     Writing graphical output: output-hill_diffusion_field_boundary_conditions/solution/solution-00002
+
+Termination requested by criterion: end time
+
+
++----------------------------------------------+------------+------------+
++----------------------------------+-----------+------------+------------+
++----------------------------------+-----------+------------+------------+
+
+-----------------------------------------------------------------------------
+-----------------------------------------------------------------------------