Make melt models work with free surface #1765

jdannberg · 2017-05-23T17:03:55Z

This is work in progress, and I basically opened the pull request so that @anne-glerum and @naliboff could have a look.
The idea is that if you have a model with melt, a new function (apply_free_surface_stabilization_with_melt) is called in the stabilization of the free surface, and it will loop over the correct DoFs (which are different for with and without melt). This is the part that I think already works.
The second point is that we also want to automatically prescribe boundary conditions for the melt velocity at the free surface, and in such a way that melt will not flow in or out of the free surface boundary. To fulfil this condition, we have to set grad p_f = rho_f g at the boundary (according to Darcy's law). We do this by not assembling the melt boundary terms at boundaries that are at the free surface (see code in melt.cc in local_assemble_stokes_system_melt_boundary).
However, we tested this and melt still seems to flow out at the boundaries, so we're currently trying to figure out why this is. But without melt being at the surface, this should already work in models with melt and free surface (the test we have is without melt at the surface).

ian-r-rose · 2017-05-23T18:05:56Z

I don't know the ins and outs of the melt system, so forgive a couple of basic questions:

You say that you are trying to set the melt velocity at the surface to have no flow in or out. Then you go on to try to set the pressure boundary condition to try to satisfy that. Is there a reason you are not directly applying velocity boundary conditions? Of course, applying both should make the system overdetermined.
What is the actual velocity boundary condition you are trying to enforce? I can think of at least three:

$v_f \cdot n = v_s \cdot n$
$v_f = v_s$
$v_f \cdot n = 0$
where s and f subscripts mean solid and fluid, respectively. Based on your description it sounds like the first one is what you are trying to implement, but I am still unsure.

jdannberg · 2017-05-23T18:13:00Z

You say that you are trying to set the melt velocity at the surface to have no flow in or out. Then you go on to try to set the pressure boundary condition to try to satisfy that. Is there a reason you are not directly applying velocity boundary conditions? Of course, applying both should make the system overdetermined.

We do not solve for the melt velocities, they are just computed as a postprocessor, using the solid velocity and the fluid pressure. So we cannot set boundary conditions for them directly.

What is the actual velocity boundary condition you are trying to enforce?

$v_f \cdot n = v_s \cdot n$

ian-r-rose · 2017-05-23T18:23:42Z

Is the fluid pressure fully determined, or is there a free constant?

ian-r-rose · 2017-05-23T18:30:51Z

I am looking at your melt paper, equation 9. It seems to me that if you set $\grad p_f = \rho_f g$, then you get $u_s = u_f$ (that is to say, option two from above).

jdannberg · 2017-05-23T18:32:48Z

There is a free constant, same as for the solid pressure. Because we solve for the fluid pressure instead of the solid pressure, we just reuse the pressure normalization we use for models without melt for the solid pressure, but apply it to the fluid pressure instead.

jdannberg · 2017-05-23T18:43:12Z

You are right, I did not write that down correctly. What we want to use as a boundary condition is that the boundary integral of some constant times $(\nabla p_f - \rho_f g) \cdot \vec{n}$ is zero.
(see equation 23 in our paper, the boundary integral on the RHS, f2 is what we can prescribe as boundary condition for the fluid pressure gradient)

ian-r-rose · 2017-05-23T18:55:29Z

Yeah, I can see how that might collide with the stabilization term. You could check this by making the stabilization parameter $\theta = zero$, thereby turning off stabilization. Timestepping would become a problem, but the velocities should be right on the first timestep. If there is still flow through the boundary, I think we could then rule out the stabilization as the cause of this particular problem.

jdannberg · 2017-05-23T20:05:37Z

That is interesting to know, thanks!
At the moment though, weirdly, what we do doesn't even work when we run the model without a free surface, and free slip or no slip boundary conditions instead (there is still a small in- or outflow).
I don't really have a good idea what's going on there...

jdannberg · 2017-06-13T16:51:42Z

@anne-glerum
In addition to the test case that is part of this pull request, I also looked at the solitary wave benchmark (in the benchmarks folder). That is in principle a 1D case, so I thought it would be easy to see if the boundary condition does what it is supposed to do. The original version of this benchmark uses the solid density times gravity as boundary fluid pressure (which is the default), so melt can flow in and out. But if you add

subsection Boundary fluid pressure model
  set Plugin name = density
  subsection Density
    set Density formulation = fluid density
  end
end

to the input file, the the solid velocity and melt velocity should be the same at the boundaries in the output (and I do not think they are at the moment, so maybe that would be a first good thing to look at).

I also added another test case to this pull request (free_surface_blob_nonzero_melt), which is principle the same as the free surface blob, just that it has a nonzero melt (so that one can check if melt flows in or out, and if the mass of porosity stays the same). I also made a simpler version of that (free_surface_blob_nonzero_melt_simple), which does not have a free surface any more, and added that input file too. In this case it is clearly visible that melt flows out at the top although it shouldn't, and if you visualize the fluid pressure gradient in paraview (Filters-->GradientOfUnstructuredDataSet), it is not equal to the product of fluid density and gravity (which is what it should be). I believe those were all of the test cases I tried.

It might also make sense to see if what I saw is just a resolution issue (if the velocity converges to the correct solution for mesh size --> 0). I remember that changing the resolution made a difference for the solution, but I did not yet do a systematic test.

What confuses me about this is that we used the fluid pressure boundary conditions for all the test cases we have in our paper (or for example the melt_transport_adaptive test in the tests directory), and we converged to the analytical solution in those cases. So before seeing the results of these tests here, I was pretty sure that our fluid pressure boundary conditions worked as intended.

jdannberg · 2018-01-26T07:53:30Z

Okay, coming back to this after a very long time...
I ran the simple version of the test (without a free surface, free_surface_blob_nonzero_melt_simple) for different resolutions, and it seems that we actually converge towards the correct solution (of a zero melt velocity at the surface), even though the convergence rate seems to be quite low. I post my results below. So I could imagine that because of the way we prescribe the boundary condition for the melt velocity (as part of the right-hand side instead of using constraints) causes the issue I saw originally: If the boundary layer is not resolved, we also do not get the correct melt velocity at the boundary. So I will continue looking at the more complicated, free surface test cases.

Refinement 3: u_f (surface) = 0.281
Refinement 4: u_f (surface) = 0.209
Refinement 5: u_f (surface) = 0.143
Refinement 6: u_f (surface) = 0.096
Refinement 7: u_f (surface) = 0.073

jdannberg · 2018-01-26T09:30:38Z

I also redid the free_surface_blob_nonzero_melt test (that has an actual free surface), setting the free surface stabilization theta to 0, and the end time to 2e4. If the velocity boundary condition for the melt works correctly, and the melt velocity is zero at the surface, melt should not flow in or out, so we can compare the global mass of the porosity after a time of 2e4 with the one at the start of the model. I post my results below, it looks like this is also converging (this time, linearly). So it seems like this is all working correctly.

Initial mass: 5.94762824e+08

Mass after 2e4 years:
5 refinements: 5.94818155e+08, difference: 5.5331e4
6 refinements: 5.94791382e+08, difference: 2.8558e4
7 refinements: 5.94775558e+08, difference: 1.2734e4
8 refinements: 5.94768854e+08, difference: 6.030e3

bangerth · 2018-01-26T14:17:21Z

So are you saying that this is now working as expected and should be merged as-is? If so, can you rebase and force push again so that the tester will run again?

jdannberg · 2018-01-27T09:51:26Z

I've done all of the testing with a rebased version of the branch that @gassmoeller made, and it seems to me that this is now working as intended. But I want to look through all the changes again to make sure everything looks good after the rebase before I force push again.

jdannberg · 2018-01-31T06:02:29Z

Okay, I pushed the rebased version and updated all of the test results, and I think everything works as intended now (see below). So it would be great if someone could have a look at this!

If melt transport is switched on, but there is no melt, we get the same results as without melt transport (only higher linear iterations numbers, which will then hopefully be fixed by our new melt solver).
If melt transport is on, and there is melt present, and we test if the global mass of the melt stays the same over time, the change in mass converges to zero with increasing resolution at the free surface (also in case of having a nonzero stabilization term for the free surface). This is both covered in the added test cases.

anne-glerum

Happy to hear that it works! I just noticed a few tiny things, otherwise it looks good to go to me!

anne-glerum · 2018-02-02T14:16:17Z

source/simulator/melt.cc

+                const Tensor<1,dim> n_hat = scratch.face_finite_element_values.normal_vector(q_point);
+                const Tensor<1,dim> g_hat = (g_norm == 0.0 ? Tensor<1,dim>() : gravity/g_norm);
+
+                double pressure_perturbation = scratch.face_material_model_outputs.densities[q_point] *


anne-glerum · 2018-02-02T14:16:27Z

source/simulator/melt.cc

+
+                const Tensor<1,dim>
+                gravity = this->get_gravity_model().gravity_vector(scratch.face_finite_element_values.quadrature_point(q_point));
+                double g_norm = gravity.norm();


anne-glerum · 2018-02-02T14:17:06Z

tests/free_surface_blob_melt.cc

+
+          for (unsigned int i=0; i < in.position.size(); ++i)
+            {
+              double porosity = std::max(in.composition[i][porosity_idx],0.0);


Make porosity const

anne-glerum · 2018-02-02T14:18:44Z

tests/free_surface_blob_nonzero_melt.prm

+# with melt transport = on and nonzero porosity.
+# The test is the same as free_surface_blob_melt, except that melt
+# has a nonzero initial condition. 
+# We test if the boundary coditions for melt are set correctly


coditions -> conditions

anne-glerum · 2018-02-02T14:20:07Z

tests/free_surface_blob_nonzero_melt.prm

+  end
+end
+
+#subsection Termination criteria


Remove lines that are commented out?

anne-glerum · 2018-02-02T14:20:44Z

tests/free_surface_blob_nonzero_melt_simple.prm

+# with melt transport = on and nonzero porosity.
+# The test is the same as free_surface_blob_melt, except that melt
+# has a nonzero initial condition. 
+# We test if the boundary coditions for melt are set correctly


coditions -> conditions

anne-glerum · 2018-02-02T14:22:46Z

tests/free_surface_blob_nonzero_melt_simple.prm

+# has a nonzero initial condition. 
+# We test if the boundary coditions for melt are set correctly
+# (no melt flowing in/out at the free surface) by checking if the
+# mass of the porosity field is conserved.


This test has no free surface (top is free slip), so the description and name are a bit confusing to me. Also there is no test output.

jdannberg · 2018-02-05T23:56:08Z

source/simulator/melt.cc

-      if (this->get_parameters().free_surface_boundary_indicators.find(boundary_indicator)
-          != this->get_parameters().free_surface_boundary_indicators.end())
-        return;
-


Just to explain this change: Even if the boundary is a free surface, someone may still want melt to flow in or out, in particular once we have the option to combine other types of surface evolution with a free surface as suggested in #2071. One example (where we noticed this would have been a problem) is modelling compaction of a partially molten layer with sedimentation from above: In this case the added pressure can squish melt out of the model.

jdannberg · 2018-02-06T00:46:57Z

Thanks for the review @anne-glerum, I think I addressed all of your comments.
So is this good to go?

jdannberg · 2018-02-09T16:45:35Z

So is this ready to merge?

gassmoeller · 2018-02-12T16:31:28Z

Looks good to me.

jdannberg force-pushed the melt_free_surface branch from ba171cc to 9702445 Compare May 23, 2017 18:22

jdannberg and others added 3 commits January 26, 2018 15:53

make free surface work with melt and the iterated IMPES scheme

4cf7282

more test cases

8200d87

Fix

f87012a

jdannberg force-pushed the melt_free_surface branch 2 times, most recently from 13ece20 to e57a6bd Compare January 30, 2018 09:17

jdannberg changed the title ~~[WIP] Make melt models work with free surface~~ Make melt models work with free surface Jan 30, 2018

more fixes and update of changelog entry

0a570b4

jdannberg force-pushed the melt_free_surface branch from e57a6bd to 0a570b4 Compare January 31, 2018 01:38

anne-glerum reviewed Feb 2, 2018

View reviewed changes

address comments

3ff1860

jdannberg commented Feb 5, 2018

View reviewed changes

gassmoeller merged commit 198ac85 into geodynamics:master Feb 12, 2018

jdannberg deleted the melt_free_surface branch February 12, 2018 21:51

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Make melt models work with free surface #1765

Make melt models work with free surface #1765

jdannberg commented May 23, 2017

ian-r-rose commented May 23, 2017

jdannberg commented May 23, 2017

ian-r-rose commented May 23, 2017

ian-r-rose commented May 23, 2017

jdannberg commented May 23, 2017

jdannberg commented May 23, 2017

ian-r-rose commented May 23, 2017

jdannberg commented May 23, 2017

jdannberg commented Jun 13, 2017

jdannberg commented Jan 26, 2018

jdannberg commented Jan 26, 2018

bangerth commented Jan 26, 2018

jdannberg commented Jan 27, 2018

jdannberg commented Jan 31, 2018

anne-glerum left a comment

anne-glerum Feb 2, 2018

anne-glerum Feb 2, 2018

anne-glerum Feb 2, 2018

anne-glerum Feb 2, 2018

anne-glerum Feb 2, 2018

anne-glerum Feb 2, 2018

anne-glerum Feb 2, 2018

jdannberg Feb 5, 2018

jdannberg commented Feb 6, 2018

jdannberg commented Feb 9, 2018

gassmoeller commented Feb 12, 2018

Make melt models work with free surface #1765

Make melt models work with free surface #1765

Conversation

jdannberg commented May 23, 2017

ian-r-rose commented May 23, 2017

jdannberg commented May 23, 2017

ian-r-rose commented May 23, 2017

ian-r-rose commented May 23, 2017

jdannberg commented May 23, 2017

jdannberg commented May 23, 2017

ian-r-rose commented May 23, 2017

jdannberg commented May 23, 2017

jdannberg commented Jun 13, 2017

jdannberg commented Jan 26, 2018

jdannberg commented Jan 26, 2018

bangerth commented Jan 26, 2018

jdannberg commented Jan 27, 2018

jdannberg commented Jan 31, 2018

anne-glerum left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jdannberg commented Feb 6, 2018

jdannberg commented Feb 9, 2018

gassmoeller commented Feb 12, 2018