"Allow fixed temperature on outflow boundaries = false" option does not work in time steps where mesh is refined

The option that allows it to prescribe the temperature only on those parts of the boundary where material flows in does not seem to work in time steps where the mesh is refined. In those time steps, the temperature is simply prescribed everywhere at the boundary. Below is an example input that reproduces the problem, where flow is prescribed to be vertically upwards and the boundary temperature in the left half of the model is slightly hotter than the initial condition.

The mesh is refined in every second step, and the boundary condition at the top boundary is also enforced every second step (even though the mesh does not actually change). The resulting problem is visible in the screenshots below:

step0 Step 0 is wrong.

step1 Step 1 is fine.

step2 Step 2 is fine.

step3 Step 3, when the mesh is refined, is wrong again.

set Dimension                              = 2

set Use years in output instead of seconds = true
set End time                               = 2e6
set Output directory                       = test

# Then there are variables that describe how the pressure should
# be normalized. Here, we choose a zero average pressure
# at the surface of the domain (for the current geometry, the
# surface is defined as the top boundary).
set Pressure normalization                 = surface
set Surface pressure                       = 0
set Adiabatic surface temperature          = 1600

subsection Geometry model
  set Model name = box
  subsection Box
    set X extent = 800e3
    set Y extent = 800e3
  end
end

subsection Initial temperature model
  set Model name = function

  subsection Function
    set Variable names      = x,y
    set Function expression = 1500 + 250 * 0.5 * (1-tanh((x-450000)/40000))
  end
end

# This is 10 K hotter than the initial condition within the plume
subsection Boundary temperature model
  set Fixed temperature boundary indicators = bottom, top
  set Allow fixed temperature on outflow boundaries = false
  set Model name = function

  subsection Function
    set Variable names      = x,y
    set Function expression = 1500 + 260 * 0.5 * (1-tanh((x-450000)/40000))
  end
end

subsection Boundary velocity model
  set Tangential velocity boundary indicators = left, right
  set Prescribed velocity boundary indicators = top: function, bottom: function

  subsection Function
    set Variable names      = x,z
    set Function expression = 0; 0.1
  end
end

subsection Gravity model
  set Model name = vertical
  subsection Vertical
    set Magnitude = 10.0
  end
end

subsection Material model
  set Model name = latent heat
  set Material averaging = project to Q1 only viscosity

  subsection Latent heat
    set Minimum viscosity             = 1e17
    set Viscosity                     = 7e21
    set Reference density             = 3291.5
    set Thermal conductivity          = 4
    set Thermal expansion coefficient = 2.6e-5
    set Reference specific heat       = 1000
    set Reference temperature         = 1600
    set Compressibility               = 0
    set Thermal viscosity exponent    = 50.0

    set Phase transition density jumps                 = 110 
    set Corresponding phase for density jump           = 0
    set Phase transition depths                        = 300000
    set Phase transition widths                        = 0
    set Phase transition temperatures                  = 1600
    set Phase transition Clapeyron slopes              = 0
    set Viscosity prefactors                           = 1,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

# The following section deals with the discretization of
# this problem, namely the kind of mesh we want to compute
# on. 
subsection Mesh refinement
  set Initial global refinement                = 4
  set Initial adaptive refinement              = 0
  set Strategy                                 = minimum refinement function
  set Time steps between mesh refinement       = 2

  subsection Minimum refinement function
    set Coordinate system = cartesian
    set Function expression = 4
  end
end

subsection Postprocess
  set List of postprocessors = visualization

  subsection Visualization
    set Time between graphical output = 0 
    set Output format                 = vtu
    set List of output variables      = material properties
  end
end

The only thing I can think of is that we forget to call compute_current_constraints() after refining. Can you check if that happens?

Thanks for helping me figure this out! I don't think that's the problem. I added print statements at the start of the refine_mesh and the compute_current_constraints() functions, and compute_current_constraints() is called in each time step (in start_timestep) before the temperature equation is solved. See output below.

-----------------------------------------------------------------------------
-- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
--     . version 2.3.0-pre (steinberger_lateral_viscosity, 21ea4ee94)
--     . using deal.II 9.4.0
--     .       with 32 bit indices and vectorization level 2 (256 bits)
--     . using Trilinos 12.18.1
--     . using p4est 2.3.2
--     . running in DEBUG mode
--     . running with 1 MPI process
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-- For information on how to cite ASPECT, see:
--   https://aspect.geodynamics.org/citing.html?ver=2.3.0-pre&sha=21ea4ee94&src=code
-----------------------------------------------------------------------------
Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 3,556 (2,178+289+1,089)

   --- Computing current constraints! --- 
*** Timestep 0:  t=0 years, dt=0 years
   --- Computing current constraints! --- 
   Solving temperature system... 0 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 30+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00000

*** Timestep 1:  t=243423 years, dt=243423 years
   --- Computing current constraints! --- 
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 15+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00001

*** Timestep 2:  t=486142 years, dt=242719 years
   --- Computing current constraints! --- 
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 15+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00002

   --- Refine mesh! --- 
Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 3,556 (2,178+289+1,089)

*** Timestep 3:  t=728580 years, dt=242439 years
   --- Computing current constraints! --- 
   Solving temperature system... 16 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 15+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00003

*** Timestep 4:  t=970878 years, dt=242298 years
   --- Computing current constraints! --- 
   Solving temperature system... 16 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 14+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00004

   --- Refine mesh! --- 
Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 3,556 (2,178+289+1,089)

*** Timestep 5:  t=1.21305e+06 years, dt=242175 years
   --- Computing current constraints! --- 
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 14+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00005

*** Timestep 6:  t=1.45511e+06 years, dt=242059 years
   --- Computing current constraints! --- 
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 13+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00006

   --- Refine mesh! --- 
Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 3,556 (2,178+289+1,089)

*** Timestep 7:  t=1.69704e+06 years, dt=241926 years
   --- Computing current constraints! --- 
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 14+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00007

*** Timestep 8:  t=1.93882e+06 years, dt=241787 years
   --- Computing current constraints! --- 
   Solving temperature system... 15 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 13+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00008

   --- Refine mesh! --- 
Number of active cells: 256 (on 5 levels)
Number of degrees of freedom: 3,556 (2,178+289+1,089)

*** Timestep 9:  t=2e+06 years, dt=61175.3 years
   --- Computing current constraints! --- 
   Solving temperature system... 11 iterations.
   Rebuilding Stokes preconditioner...
   Solving Stokes system... 11+0 iterations.

   Postprocessing:
     Writing graphical output: test/solution/solution-00009

Termination requested by criterion: end time

+----------------------------------------------+------------+------------+
| Total wallclock time elapsed since start     |      13.5s |            |
|                                              |            |            |
| Section                          | no. calls |  wall time | % of total |
+----------------------------------+-----------+------------+------------+
| Assemble Stokes system           |        10 |      3.01s |        22% |
| Assemble temperature system      |        10 |      4.19s |        31% |
| Build Stokes preconditioner      |        10 |      2.25s |        17% |
| Build temperature preconditioner |        10 |    0.0152s |      0.11% |
| Initialization                   |         1 |     0.678s |         5% |
| Postprocessing                   |        10 |     0.896s |       6.7% |
| Refine mesh structure, part 1    |         4 |      0.29s |       2.2% |
| Refine mesh structure, part 2    |         4 |    0.0389s |      0.29% |
| Setup dof systems                |         5 |     0.139s |         1% |
| Setup initial conditions         |         1 |     0.129s |      0.96% |
| Setup matrices                   |         9 |     0.602s |       4.5% |
| Solve Stokes system              |        10 |     0.142s |       1.1% |
| Solve temperature system         |        10 |    0.0159s |      0.12% |
+----------------------------------+-----------+------------+------------+

It must be something about the velocity values being stored in the current linearization point. When I replace current_linearization_point by old_solution in the replace_outflow_boundary_ids function here https://github.com/geodynamics/aspect/blob/07252b14f88d763cdcf4db6cb02238fcd97ebfe6/source/simulator/helper_functions.cc#L2105 then it works.

Are we forgetting to apply current_constraints with the correct velocity boundary conditions to the solution vector?

I think the problem is that we call compute_current_constraints() before we call initialize_current_linearization_point(). Or at least I tested that if I add a line with initialize_current_linearization_point() before compute_current_constraints(), then the temperature boundary conditions are applied correctly. Just to demonstrate what I mean, I have opened a PR here: #5043.

So it seems that something is wrong with the values being stored in the current_linearization_point if we refine the mesh without calling initialize_current_linearization_point(). Does my suggestion make sense or is there a better way of fixing this?

In addition, the temperature boundary conditions are always applied in time step 0 (I think because we have not solved for the velocity yet, and a zero velocity does not count as outflow; we check that the integrated outwards flow through a boundary face > 0, so if the flow is zero, that's not counted as outflow). This also does not quite make sense (at least if we know the boundary velocity, we should be able to use that rather than a zero velocity) and it would be great to find a fix for that as well.

geodynamics / aspect

"Allow fixed temperature on outflow boundaries = false" option does not work in time steps where mesh is refined #5042