Account for the Darcy Velocity and/or the Fluid Velocity when Determining Outflow Boundaries

[danieldouglas92]

Currently, the replace_outflow_boundary_ids function only uses the solid velocity to determine if a boundary should be flagged as an outflow boundary. This PR updates the replace_outflow_boundary_ids function to account for the darcy velocity when using the darcy field compositional field advection method.

Based on a conversation with @jdannberg this will also need to be extended to account for the fluid velocity when simulating melt transport, but I figured I would open the PR with what I have first to make sure that the implementation looks good and then add the other case afterwards.

This PR addresses the last point in Issue #4748.

[danieldouglas92]

Thanks for the quick review @tjhei! I've added another if statement which accounts for the fluid velocity used by the melt transport model when determining whether a boundary is an outflow boundary

[jdannberg]

Thank you for taking a stab at this!
I think this goes in the right direction. But I think it would be good to have a test case that shows that this actually works (one with darcy flow where the solid flows in and the melt flow out, for example).

[jdannberg]

I know you said you wanted to split this into 2 PRs and this will be changed later when you actually take into account the different field types. But I think like this, it might break some models. Because now, if melt flows out, and solid flows in, you would not set boundary conditions. But you need boundary conditions when material flows in (otherwise, what temperature/composition will the inflowing material have?) So at the very least only change it to an outflow boundary if all fields flow out (which is a bit tricky with how you've set up the integrated_flow variables right now).

[danieldouglas92]

This makes sense! I've changed this so that they all need to be outflowing to update the cell.

[danieldouglas92]

Thanks @jdannberg for the review, I'll add a test case in a separate commit soon!

[tjhei]

these two lines seem redundant.

[danieldouglas92]

I do this because by default (no melt transport) integrated_fluid_flow = std::numeric_limits<double>::max() so that the if statement at line 2400 still works in solid only models. If there is melt though, then I need to set integrated_fluid_flow = 0 otherwise line 2388:

integrated_fluid_flow += (face_current_fluid_velocity_values[q] * fe_face_values.normal_vector(q)) *
                         fe_face_values.JxW(q);

will always be infinity. This was the only way I could think to get around this issue, but if you have a better approach that would be great!

[danieldouglas92]

The same is done for when a field is advected using the darcy advection method at line 2352

[tjhei]

Oh, I see. I would put that into the initialization directly then by using a ternary (?) operator.

[danieldouglas92]

Ah yes! That will make this look much cleaner

[danieldouglas92]

@jdannberg @tjhei John (@naliboff) and I just talked about what this PR is now doing given that the condition for flagging a cell as an outflow cell requires outflow of all three velocities and we're wondering if we should just scrap this PR and jump right into making it composition dependent.

The issue that this PR was trying to address was that in a model with two-phases the solid phase could be inflowing while the fluid phase is outflowing, this is why I had initially implemented line 2388 as an OR condition and not an AND condition. Since it is an AND condition now (to not violate the boundary condition for the solid as it inflows), the original issue that this PR is trying to address is no longer solved.

Additionally, in the rare case where the solid is outflowing and the fluid is inflowing, the cell will not be flagged as outflow because both phases are not outflowing, which is an issue that wasn't present prior to this PR.

Am I missing something that makes what I'm currently trying to do still valuable? If not, I will close this PR and work on a composition dependent implementation.

[jdannberg]

I think it makes sense to directly go to implementing this in dependence of the compositional field type. But I think you can just work from what you have here. The only changes that would be required are

(1) add an input parameter to the replace_outflow_boundary_ids function that tells the function if it computes boundary ids for the temperature or compositional field, and which field,
(2) to modify what you have right now so that it checks for the compositional field type (or temperature), and if it's the porosity (or a melt field) in a model with melt transport, or a darcy field, you pick the appropriate velocity to determine if model has outflow (rather than the current AND/OR statement)
(3) change core.cc so that replace_outflow_boundary_ids and restore_outflow_boundary_ids is called within the loop over all compositional fields (Do this and this in a loop for all compositional field and also for the temperature)

[danieldouglas92]

@jdannberg Here's my initial go at making this composition dependent. I wasn't sure if I needed to make restore_outflow_boundary_ids composition dependent so I left this function as it was for now to see how the testers did.

[jdannberg]

I think this should work!

There's one case (porosity in the case of melt transport) which I think does not work correctly right now (see other comments), but I think in general this should now apply the different velocities correctly at the boundary. You also do not need to change restore_outflow_boundary_ids because that just resets the values to what they were originally.

Did you already check the tests to see if the results look right? Because it would make sense that some tests are different now (since you changed the functionality).

[danieldouglas92]

Thanks for the review @jdannberg! The tests are failing because the boundary_ids are exceeding 128, working on fixing this bug now but hopefully after that is fixed this is good to go!

[danieldouglas92]

Think I fixed it! restore_outflow_boundary_ids was only being called if the composition discretization was continuous, working on a testcase now

[danieldouglas92]

While trying to come up with a test I ran into an issue with the case where consider_darcy_velocity = true. When this part of the code:

                if (consider_darcy_velocity)
                  {
                    in.reinit(fe_face_values, cell, introspection, solution);
                    material_model->evaluate(in, out);
                    fluid_out = out.template get_additional_output<MaterialModel::MeltOutputs<dim>>();
                  }

gets called (specifically material_model->evaluate(in, out);), the model terminates because the temperature, pressure etc. are all 0.

A similar error happens for the same test if porosity is not advected with the darcy field advection method, but instead with the fluid velocity calculated by the melt model. In core.cc at when this line gets called, this triggers the same error when the code this->get_material_model().evaluate(material_model_inputs, material_model_outputs); in melt.cc gets called and terminates due to the temperature, pressure etc. all being 0.

I'm not sure what to make of this, as for the fluid velocity case this all happening in functions that I didn't modify in this PR. The only connection between what I changed and the fluid velocity case is through the variable new_current_constraints, but it's not clear to me how what I've done would propagate through into this variable. I've added a commit just now which includes the test that I was working on in an unfinished state. I would appreciate any guidance on how to proceed here!

[jdannberg]

I think I don't yet understand what exactly the problem is you're describing.
What exactly is the error that makes the model crash? In other words, why is it a problem if temperature, pressure etc. are all zero? (I get that it's not what you put in, but that's a different issue.) At which point does it crash? (for example, is the temperature zero because it's before the initial conditions have been set, or should they have been set and for some reason are zero now?)

The line
in.reinit(fe_face_values, cell, introspection, solution)
just takes whatever values are in solution and puts those into in, so what do you expect to be in the solution vector at that point in time where it crashes and is that the actually the case?

Does your input file crash if you run it on ASPECT main? Or is the crash actually related to the changes you made?

[danieldouglas92]

@jdannberg Yes sorry for not including more background detail on this. My first instinct was, as you pointed out could be an issue, was that I was trying to use the material model inputs before they were initialized and so were set to 0.

For the Darcy velocity case, this error does only happens on this branch, and does not happen on main. For the fluid velocity case it happens on both main and on my branch. It isn't surprising that this doesn't happen for the Darcy velocity case on main though, because the line material_model->evaluate(in, out); isn't there.

I think I'm just now figuring out why this is happening. The test no_dirichlet_on_outflow_melt, and most of the tests involving melt, use the melting rate material model fromrising_melt_blob.cc which explicitly defines what the outputs should be regardless of the modeled temperature and pressure. This test that I've created relies on using the visco plastic material model to determine the material outputs based on the pressure and temperature, but these are not initialized yet when replace_outflow_boundary_ids is called.

Of course I could probably circumnavigate this by using rising_melt_blob.cc (I'll double check that this actually works soon) in this test, but this wouldn't actually fix the issue. If you agree that my view on the problem is right, what do you think would be the best way to get around this? I have some ideas for using a reference pressure or temperature only when the function is called per-initialization of the material model inputs, but they seem kind of hacky.

[danieldouglas92]

@jdannberg I got around the issue by properly setting the initial adiabatic conditions! Now the tests are not crashing, but they are producing some behaviour that is really bizarre. For both cases, the test simulates a blob of porosity in the middle of the model which flows up to outflow through the top boundary. At the top boundary, there is a solid field called solid_phase which is being advected downward at 10 cm/yr due to the prescribed velocity boundary conditions. The hope is that this will allow the solid_phase composition to be correctly handled as an inflow while the porosity composition will be allowed to outflow at the same time. This is what the initial conditions look like

For the fluid velocity case, when the blob hits the top boundary the porosity starts to increase across the whole domain and the model crashes. This is what happens to the porosity composition:

For the Darcy velocity case though, this seems to be ok.

I'm not exactly sure what is causing the fluid velocity case to crash. The model hangs when run in debug mode, and terminates with an MPI error in release mode.

[jdannberg]

For the first problem with the pressure and temperature: Are you saying you need to set the adiabatic conditions to a function, and if you use initial profile, it does not work? If that’s the case, we have many places in the material models where we check if this->get_adiabatic_conditions().is_initialized() and then do something different if it is not initialized. So whichever computation of material properties can’t deal with adiabatic pressure and temperature being unavailable, that’s what would have to be changed. But that does not directly relate to your PR, it’s just something we should fix independently. Since you say it’s broken on main, could you make an issue describing the problem in detail?

For the blob: That does seem strange. And I would say, even for the Darcy velocity, it does not seem right. It looks like there is a small blue line at the top of the box, so the outflow boundary conditions are not actually recognized correctly and porosity is still set to zero at the boundary (whereas at the bottom boundary, porosity along the boundary is the same as what it is in the interior).

For the fluid velocity case: where does it hang in debug mode? Can you run it in a debugger and figure out the line? And it looks like you’ve set up the model to not have melting or freezing, so it really just advects the porosity field? Is there anything else that is weird about the model (pressure, temperature, etc.)? Does it also look like that and crash on main or is it related to the boundary conditions?

And at least the boundary condition part might work for the porosity? It is set to zero along parts of the upper and the whole lower boundary, but not where the blob reached the boundary.

[danieldouglas92]

@jdannberg Thanks for following up on this!

For the first problem, yes even on the main branch if I do not define an Adiabatic conditions model (which then uses the compute profile model as the default), the model immediately crashes because the temperature and pressure are both 0. When I set the Adiabatic conditions model with a function which explicitly defines the pressure and temperature, then the model proceeds. (PS I realized that this is because I wasn't explicitly defining an Adiabatic surface temperature, and I had just assumed the default value would be nonzero, but the default value is in fact 0. Setting the adiabatic surface temperature overcomes this error).

For the darcy flow case, you are right I didn't notice that there is a thin blue line along the top boundary where the porosity is trying to outflow. I wonder if the bottom boundary looking correct might not be surprising since the solid is also outflowing here?

For the fluid velocity case, I wasn't able to found out an exact line unfortunately, but I was
able to narrow it down to this when this function is called here. When I run it with the debugger it makes it to the end of this function when advecting the porosity composition, and then it throws a Trilinos exception. This is the exact message:

---------------------------------------------------------
TimerOutput objects finalize timed values printed to the
screen by communicating over MPI in their destructors.
Since an exception is currently uncaught, this
synchronization (and subsequent output) will be skipped
to avoid a possible deadlock.
---------------------------------------------------------
2

----------------------------------------------------
Exception 'ExcTrilinosError(ierr)' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <1678> of file </Users/danieldouglas/software/dealii-9.5.2-install/tmp/unpack/deal.II-v9.5.2/source/lac/trilinos_sparse_matrix.cc> in function
    void dealii::TrilinosWrappers::SparseMatrix::add(const size_type, const size_type, const size_type *, const TrilinosScalar *, const bool, const bool)
The violated condition was: 


For historical reasons, many Trilinos functions express errors by
returning specific integer values to indicate certain errors.
Unfortunately, different Trilinos functions often use the same integer
values for different kinds of errors, and in most cases it is also not
documented what each error code actually means. As a consequence, it
is often difficult to say what a particular error (in this case, the
error with integer code '2') represents and how one should fix a code
to avoid it. The best one can often do is to look up the call stack to
see which deal.II function generated the error, and which Trilinos
function the error code had originated from; then look up the Trilinos
source code of that function (for example on github) to see what code
path set that error code. Short of going through all of that, the only
other option is to guess the cause of the error from the context in
which the error appeared.

Unfortunately I haven't been able to make much of this error message... I found the error message in the dealii source code but as the error message points out I didn't find it very helpful. If you've encountered this before or have any ideas on how to proceed I would really appreciate it.

Also, here is the stacktrace in case that might be helpful:

Stacktrace:
-----------
#0  2   libdeal_II.g.9.5.2.dylib            0x0000000113efae00 _ZN6dealii16TrilinosWrappers12SparseMatrix3addEjjPKjPKdbb + 532: 2   libdeal_II.g.9.5.2.dylib            0x0000000113efae00 _ZN6dealii16TrilinosWrappers12SparseMatrix3addEjjPKjPKdbb 
#1  3   libdeal_II.g.9.5.2.dylib            0x0000000113bde68c _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_16TrilinosWrappers17BlockSparseMatrixENS3_3MPI11BlockVectorEEEvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjEEEERT_RT0_bNSF_17integral_constantIbLb1EEE + 1156: 3   libdeal_II.g.9.5.2.dylib            0x0000000113bde68c _ZNK6dealii17AffineConstraintsIdE26distribute_local_to_globalINS_16TrilinosWrappers17BlockSparseMatrixENS3_3MPI11BlockVectorEEEvRKNS_10FullMatrixIdEERKNS_6VectorIdEERKNSt3__16vectorIjNSF_9allocatorIjEEEERT_RT0_bNSF_17integral_constantIbLb1EEE 
#2  4   aspect-debug                        0x00000001004ecdd4 _ZN6aspect9SimulatorILi2EE37copy_local_to_global_advection_systemERKNS1_14AdvectionFieldERKNS_8internal8Assembly8CopyData15AdvectionSystemILi2EEE + 92: 4   aspect-debug                        0x00000001004ecdd4 _ZN6aspect9SimulatorILi2EE37copy_local_to_global_advection_systemERKNS1_14AdvectionFieldERKNS_8internal8Assembly8CopyData15AdvectionSystemILi2EEE 
#3  5   aspect-debug                        0x0000000100500c2c _ZN6dealii10WorkStream8internal10sequential3runIZN6aspect9SimulatorILi2EE25assemble_advection_systemERKNS6_14AdvectionFieldEEUlRKNS_18TriaActiveIteratorINS_15DoFCellAccessorILi2ELi2ELb0EEEEERNS4_8internal8Assembly7Scratch15AdvectionSystemILi2EEERNSH_8CopyData15AdvectionSystemILi2EEEE_ZNS6_25assemble_advection_systemES9_EUlRKSO_E_NS_16FilteredIteratorISD_EESK_SO_EEvRKT1_RKNS_9std_cxx2013type_identityISW_E4typeET_T0_RKT2_RKT3_ + 212: 5   aspect-debug                        0x0000000100500c2c _ZN6dealii10WorkStream8internal10sequential3runIZN6aspect9SimulatorILi2EE25assemble_advection_systemERKNS6_14AdvectionFieldEEUlRKNS_18TriaActiveIteratorINS_15DoFCellAccessorILi2ELi2ELb0EEEEERNS4_8internal8Assembly7Scratch15AdvectionSystemILi2EEERNSH_8CopyData15AdvectionSystemILi2EEEE_ZNS6_25assemble_advection_systemES9_EUlRKSO_E_NS_16FilteredIteratorISD_EESK_SO_EEvRKT1_RKNS_9std_cxx2013type_identityISW_E4typeET_T0_RKT2_RKT3_ 
#4  6   aspect-debug                        0x00000001004edc24 _ZN6dealii10WorkStream3runIZN6aspect9SimulatorILi2EE25assemble_advection_systemERKNS4_14AdvectionFieldEEUlRKNS_18TriaActiveIteratorINS_15DoFCellAccessorILi2ELi2ELb0EEEEERNS2_8internal8Assembly7Scratch15AdvectionSystemILi2EEERNSF_8CopyData15AdvectionSystemILi2EEEE_ZNS4_25assemble_advection_systemES7_EUlRKSM_E_NS_16FilteredIteratorISB_EESI_SM_EEvRKT1_RKNS_9std_cxx2013type_identityISU_E4typeET_T0_RKT2_RKT3_jj + 220: 6   aspect-debug                        0x00000001004edc24 _ZN6dealii10WorkStream3runIZN6aspect9SimulatorILi2EE25assemble_advection_systemERKNS4_14AdvectionFieldEEUlRKNS_18TriaActiveIteratorINS_15DoFCellAccessorILi2ELi2ELb0EEEEERNS2_8internal8Assembly7Scratch15AdvectionSystemILi2EEERNSF_8CopyData15AdvectionSystemILi2EEEE_ZNS4_25assemble_advection_systemES7_EUlRKSM_E_NS_16FilteredIteratorISB_EESI_SM_EEvRKT1_RKNS_9std_cxx2013type_identityISU_E4typeET_T0_RKT2_RKT3_jj 
#5  7   aspect-debug                        0x00000001004ed634 _ZN6aspect9SimulatorILi2EE25assemble_advection_systemERKNS1_14AdvectionFieldE + 1232: 7   aspect-debug                        0x00000001004ed634 _ZN6aspect9SimulatorILi2EE25assemble_advection_systemERKNS1_14AdvectionFieldE 
#6  8   aspect-debug                        0x0000000100103b34 _ZN6aspect9SimulatorILi2EE30assemble_and_solve_compositionERKNSt3__16vectorIdNS2_9allocatorIdEEEEPS6_ + 936: 8   aspect-debug                        0x0000000100103b34 _ZN6aspect9SimulatorILi2EE30assemble_and_solve_compositionERKNSt3__16vectorIdNS2_9allocatorIdEEEEPS6_ 
#7  9   aspect-debug                        0x0000000100105188 _ZN6aspect9SimulatorILi2EE35solve_iterated_advection_and_stokesEv + 336: 9   aspect-debug                        0x0000000100105188 _ZN6aspect9SimulatorILi2EE35solve_iterated_advection_and_stokesEv 
#8  10  aspect-debug                        0x0000000100ae4998 _ZN6aspect9SimulatorILi2EE14solve_timestepEv + 252: 10  aspect-debug                        0x0000000100ae4998 _ZN6aspect9SimulatorILi2EE14solve_timestepEv 
#9  11  aspect-debug                        0x0000000100ae3498 _ZN6aspect9SimulatorILi2EE3runEv + 892: 11  aspect-debug                        0x0000000100ae3498 _ZN6aspect9SimulatorILi2EE3runEv 
#10  12  aspect-debug                        0x0000000100a633e8 _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbbb + 1076: 12  aspect-debug                        0x0000000100a633e8 _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbbb 
#11  13  aspect-debug                        0x0000000100a62a40 main + 2300: 13  aspect-debug                        0x0000000100a62a40 main 
#12  14  dyld                                0x000000018ea0b154 start + 2476: 14  dyld                                0x000000018ea0b154 start 
--------------------------------------------------------

[jdannberg]

Oh, I see. So everything works fine with the adiabatic conditions as long as you set the adiabatic surface temperature. That's good to hear.

For the fluid velocity case, for me it worked with a high resolution everywhere:

subsection Mesh refinement
  set Initial adaptive refinement              = 0
  set Initial global refinement                = 2
  set Time steps between mesh refinement       = 0
end

I am still not sure what the problem is, and it would be good to figure that out, but you could try if your current file also crashes on main, and if so, then at least you know it does not have anything to do with your implementation. Since I had simplified a few more things, I am attaching my input file here.
original.prm.txt

I also still see the porosity increasing in the box at the end even though I don't think any fluid flows in and there are no reactions. That still seems quite wrong, but at the moment I can't see where that comes from.

Looking at the output, for the fluid flow case the boundary conditions seem correct for the porosity (for the solid_phase field it is hard to tell, because the value at the outflow boundary is the same as the initial condition).

I wonder if the bottom boundary looking correct might not be surprising since the solid is also outflowing here?

Yes, I think that makes sense and gives you an idea where to search for the bug.

[danieldouglas92]

@jdannberg I'm pretty confident I've figured out what is going on with these tests. First, when I increase the resolution of the fluid velocity test, the behaviour does drastically improve, and the behaviour of the test using the functionality implemented in this PR relative to the main branch definitely suggests that it is working as expected.

On the main branch, because the top boundary is strictly an inflow boundary, the porosity pools up and spreads out along the top boundary. With the functionality included in this PR, the porosity is allowed to outflow while the solid inflows, and no pooling happens at the top boundary. Two screenshots from the same model time comparing the behaviour on main vs this PR are shown below.

For the Darcy velocity test, the issue is that the Darcy velocity depends on the porosity compositional field. In this test, when I had set that the composition should be fixed on the top and bottom boundaries to the initial composition, the initial composition has the porosity field equal to 0 on both of these boundaries. This means that the darcy velocity will ALWAYS be essentially whatever the solid velocity is, and in this test case since the solid velocity is NEGATIVE on the top boundary, the porosity will NEVER outflow. The only way to allow the porosity to outflow is for there to BE porosity, and it's never able to become non-zero due to the enforced boundary conditions.

If I define the porosity to be uniform everywhere and create an additional field which is advected with the darcy field method (following the setup of tests/darcy_velocity.prm), the model behaves as expected and the behaviour is more correct relative to the main branch. Below are a few of screenshots comparing the behaviour of this test between the changes in this PR and with main.

The first pair of screenshots shows how when using the main branch the free_fluid composition is zero at the top boundary since the boundary is strictly an inflow boundary, and also how numerical instabilities arise later in the model.
Note the colorbar scale change in the second screenshot, where the composition blows up.

The second pair of screenshots shows the behaviour on this branch, where these issues are not present. The composition effectively outflows through the top boundary and avoids blowing up later in the model.

While I'm confident that this is why the Darcy velocity test was not working, the solution of defining a homogeneous porosity is obviously not ideal. The framework that we've put together for using the darcy field advection method with the Reactive Fluid Transport Model requires the porosity compositional field to be dynamic and most models likely will not start with non-zero porosity on all boundaries. I'm not quite sure how to get around this issue though, so I would love to hear your thoughts on this Juliane!