Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Bugfix] add an option that fixes temperature/composition only on inflow/close… #2419

Merged
merged 13 commits into from Oct 6, 2018

Conversation

jdannberg
Copy link
Contributor

…d boundaries, but not on outflow boundaries.

Fixes #903.

We can discuss what the default for the temperature equation should be.

And I still have to add a changelog entry (which I will do after we've agreed on a the default for the temperature).

// Now check if the current_constraints we just computed changed from before. Do this before the melt handler
// adds constraints for the melt cells, because they are allowed to change (and often do) without us having
// to reinit the sparsity pattern.
bool constraints_changed = false;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

move this inside the if (!mesh_has_changed)

bool mesh_has_changed = false;
if (!(current_constraints.get_local_lines().size() == new_current_constraints.get_local_lines().size()))
mesh_has_changed = true;
else if (!(current_constraints.get_local_lines() == new_current_constraints.get_local_lines()))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if ((A.size() != B.size()) || (A != B))
  mesh_has_changed =true;

current_constraints.clear ();
current_constraints.reinit (introspection.index_sets.system_relevant_set);
current_constraints.merge (new_current_constraints);
#endif
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add current_constraints.close(); before endif

current_constraints.clear ();
current_constraints.reinit (introspection.index_sets.system_relevant_set);
current_constraints.merge (new_current_constraints);
#endif

if (time_step == 0 && parameters.include_melt_transport)
melt_handler->add_current_constraints (current_constraints);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am worried that the melt handler will break our logic: if we add constraints here, we will rebuild the sparsity pattern even though that shouldn't be necessary. I am not quite sure how to fix this. Maybe ignore p_c indices in the check above?!

@jdannberg
Copy link
Contributor Author

There were some tests failing, and I mostly modified the input files so that they now exactly reproduce the old behaviour. The only change in the tests I updated was the memory consumption of the constraints, so I guess this is expected.

I also addressed the comments. Can someone have another look?

Copy link
Contributor

@bangerth bangerth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I read through the first part, but it's late. I'll check the rest tomorrow, after you've already addressed some of the initial comments.

* on parts of the boundaries where material flows out.
*/
bool
allows_fixed_composition_on_outflow_boundaries() const;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All other names here are in the infinitive form, like in find_boundary_... or get_fixed_.... Maybe change to allow_fixed_...?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I modeled this after other functions we have that return bools, like has_boundary_temperature, simulator_is_initialized, pressure_rhs_needs_compatibility_modification etc. If the function returns a bool, I think this naming scheme makes sense, as you would use it in statements like
if (boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries()).
That seems more intuitive than
if (boundary_composition_manager.allow_fixed_composition_on_outflow_boundaries()).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see -- it's a predicate (you're asking a true/false question) rather than a request for information. OK then.

* Whether we allow the composition to be fixed on parts of the boundary
* where material flows out of the domain.
*/
bool allow_fixed_composition_on_outflow_boundaries;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You'll have to change this name then since it'll clash with the changed function name.

* on parts of the boundaries where material flows out.
*/
bool
allows_fixed_temperature_on_outflow_boundaries() const;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same here, obviously

/**
* Offset the boundary id of all faces located on an outflow bounday
* by a fixed value given by the input parameter
* @param boundary_id_offset.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dislike that this is a parameter. That's something that will be completely bizarre to basically every user who does not look up the code.

I'd much rather than you simply disallow a certain range of boundary indicators (and add a check somewhere) and then do the shift internally.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is how I had it in the beginning, but I need the offset in two functions (replace_outflow_boundary_ids and reset_outflow_boundary_ids), and it has to be the same in both. This way I at least make sure that the same number is used in both places. Do you have a good suggestion on how to do this?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you already hard-coded it. That's ok. I know how to fix it later on.

@@ -1193,6 +1193,25 @@ namespace aspect
LinearAlgebra::BlockVector &relevant_dst,
LinearAlgebra::BlockVector &tmp_distributed_stokes);

/**
* Offset the boundary id of all faces located on an outflow bounday
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

-> boundary

* This function is implemented in
* <code>source/simulator/helper_functions.cc</code>.
*/
void reset_outflow_boundary_ids(const unsigned int boundary_id_offset);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think "restore" would be a better word.

"When the composition is fixed on a given boundary as determined "
"by the list of 'Fixed composition boundary indicators', there "
"might be parts of the boundary where material flows out and "
"one may want to prescribe the composition only on the parts of "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the -> those

"by the list of 'Fixed composition boundary indicators', there "
"might be parts of the boundary where material flows out and "
"one may want to prescribe the composition only on the parts of "
"the boundary where there is inflow. This parameter determines "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the -> those

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this comment (but fixed the one above for line 334). Was this what you meant?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's just a question of emphasis. "only on the parts the boundary" and "only on those parts of the boundary" means the same, it just provides additional emphasis.

"the boundary where there is inflow. This parameter determines "
"if compositions are only prescribed at these inflow parts of the "
"boundary (if false) or everywhere on a given boundary, independent "
"of the flow direction (if true).");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add the following paragraph to this:

"Mathematically speaking, the compositional fields satisfy an advection equation that has no diffusion. For this equation, one can only impose Dirichlet boundary conditions (i.e., prescribe a fixed compositional field value at the boundary) at those boundaries where material flows in. This would correspond to the true'' setting of this parameter, which is correspondingly the default. On the other hand, on a finite dimensional discretization such as the one one obtains from the finite element method, it is possible to also prescribe values on outflow boundaries, even though this may make no physical sense. This would then correspond to the false'' setting of this parameter."

"the boundary where there is inflow. This parameter determines "
"if temperatures are only prescribed at these inflow parts of the "
"boundary (if false) or everywhere on a given boundary, independent "
"of the flow direction (if true).");
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

here add something similar, except that it is necessary to state that mathematically it is allowed to prescribe the temperature even on outflow boundaries as long as the diffusion coefficient is nonzero. In practice, however, this would only make physical sense if the diffusion coefficient is actually quite large to prevent the creation of a boundary layer.

@jdannberg
Copy link
Contributor Author

@bangerth, I think I addressed all of your comments, so it would be great if you could have another look.

Copy link
Contributor

@bangerth bangerth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One comment before I get to the details: I think "Allow fixed composition/temperature" is maybe not the right term for the parameter/member variable/member function because the field does not have to be fixed -- which to me would imply that it is constant in time. How about you replace "fixed" by "prescribed"?

* on parts of the boundaries where material flows out.
*/
bool
allows_fixed_composition_on_outflow_boundaries() const;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see -- it's a predicate (you're asking a true/false question) rather than a request for information. OK then.

@jdannberg
Copy link
Contributor Author

For the name of the parameter, I just wanted to make that consistent with the rest of the parameters/code. For the temperature and composition, we call the input parameter Fixed temperature boundary indicators and Fixed composition boundary indicators (and the same is true for the names of the variables in the code). You can argue about how much sense those names make, but my new parameter refers to these old names. So I think "Allow fixed composition/temperature on outflow boundaries" is the right choice here.

@bangerth
Copy link
Contributor

Ah, so much historical baggage :-( In that case, leave it as you have it but maybe make a note in the doc of the parameter that fixed = Dirichlet, not time independent.

Copy link
Contributor

@bangerth bangerth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here my remaining comments.

Let's postpone the discussion about the replace/restore functions to a follow-up. I can do this myself; I think that will be faster than trying to explain how I would like that function to look like.

// we update the boundary indicators of all faces that belong to ouflow boundaries
// so that they are not in the list of fixed composition boundary indicators any more.
// We will undo this change in a later step, after the constraints have been set.
if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries())
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a copy of the comment from above except that you want to apply it for the compositional fields now. maybe say something along the lines "Use the same trick for marking up outflow boundary conditions for compositional fields as we did above already for the temperature."

I guess you convinced yourself that you need the correct boundary indicators in the call to update above in line 751 (new) and that you can't just avoid the restore/replace pair in line 747/758?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, I just thought it would be cleaner like this. There are different parameters in the if-statement, so we would have to do something like
if (!boundary_temperature_manager.allows_fixed_temperature_on_outflow_boundaries() && boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries()) restore_outflow_boundary_ids(boundary_id_offset); if (!boundary_composition_manager.allows_fixed_composition_on_outflow_boundaries() && boundary_temperature_manager.allows_fixed_temperature_on_outflow_boundaries()) replace_outflow_boundary_ids(boundary_id_offset);

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I missed that. I thought you're just restoring it for the update call immediately above.

OK then, just edit the comment.


// let the melt handler add its constraints once before we solve the porosity system for the first time
if (parameters.include_melt_transport)
melt_handler->save_constraints (current_constraints);
melt_handler->save_constraints (new_current_constraints);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is pre-existing, but it still makes me unsure: the function is called "save_" which suggests that the melt handler stores a copy somewhere. Is this correct? Does it copy, or does it keep a reference? Keeping a reference will not work any more, of course.

I wonder about this because the comment in front of the line does not seem to be in sync with the code -- it talks about adding constraints, when the function is called "save_". There is also a comment a couple of lines below in which you say "Do this before the melt handler adds constraints...", which seems to be in conflict with the comment here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The melt handler copies the current constraints and stores them in a private member variable so that we can use them later. This is necessary because we want to add the melt constraints, which depend on the solution of the porosity field, later on, after we have computed this solution. In this way, we only need to update the constraints matrix instead of computing all constraints again.
(at least this is what the documentation in melt.h says)

But you are right, the comment refers to a different line, and I'll update it.

// If the mesh got refined and the size of the linear system changed, the old and new constraint
// matrices will have different entries, and we can not easily compare them.
// However, in case of mesh refinement we have to rebuild the matrix anyway, so we can skip the
// step that checks if constraints have changed.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

couldn't you do this easier by invalidating the current_constraints object at the time of mesh refinement?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@tjhei designed this part and can answer that question much better than I can.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's not worry about it then.

// step that checks if constraints have changed.
bool mesh_has_changed = false;
if (!(current_constraints.get_local_lines().size() == new_current_constraints.get_local_lines().size())
|| !(current_constraints.get_local_lines() == new_current_constraints.get_local_lines()))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you really need both comparisons? I would have expected that if you ask for the two objects to be the same (second line), that the sizes must also be the same (first line).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, but if the sizes are different (first line) the code will crash if I just call the second line.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right.

I think you could still make this easier to read by using a != b instead of !(a==b).

bool mesh_has_changed = false;
if (!(current_constraints.get_local_lines().size() == new_current_constraints.get_local_lines().size())
|| !(current_constraints.get_local_lines() == new_current_constraints.get_local_lines()))
mesh_has_changed = true;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can roll the condition you use in the if here right into the initialization of the mesh_has_changed variable.

#if DEAL_II_VERSION_GTE(9,0,0)
for (auto &row: current_constraints.get_lines())
{
if (!new_current_constraints.is_constrained(row.index))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All you check here is whether the set of constrained DoFs is the same, but not whether the constraint weights are the same. Is this on purpose? If so, maybe not call the variable constraints_changed but constrained_dofs_set_changed?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. Can you leave a note, @jdannberg ? This is to decide if we need to construct a new sparsity pattern. The matrices will be reassembled in each timestep regardless and the value of the constraints does not matter for the sparsity pattern.

for (const auto row: current_constraints.get_local_lines())
{
if (current_constraints.is_constrained(row)
!= new_current_constraints.is_constrained(row))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With these conditions here (for 8.5, and the corresponding one for 9.0), you really only change that the set of constrained DoFs we have currently is a subset of the ones we had before. But there may have been more before. Is that what you want to test? If you really want equality, then you also need to check the other way around.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we concluded we don't have to check the other way around because if you only add constraints you don't have to rebuild the matrix. But I may have gotten that wrong, @tjhei, can you add to that?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That seems wrong. New constraints may lead to new entries in the matrix.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree that both versions should do the same thing (even though either way will work)

"if compositions are only prescribed at these inflow parts of the "
"boundary (if false) or everywhere on a given boundary, independent "
"of the flow direction (if true)."
"Mathematically speaking, the compositional fields satisfy an "
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

add "\n\n" to provide offset here.

maybe do the same for the corresponding temperature entry.

@jdannberg
Copy link
Contributor Author

Okay, I think I addressed all of @bangerth's comments, except for the three marked for @tjhei. So @tjhei, maybe you can have a look?

@jdannberg
Copy link
Contributor Author

Now that I think about it, for the porosity and any other fields that are advected with the melt velocity, we should probably use the melt velocity to decide if the boundary is an outflow boundary, correct?

@gassmoeller
Copy link
Member

Yes, whatever velocity is used in the advection term.

@jdannberg
Copy link
Contributor Author

Okay, I thought about this, and there's another problem with that. In most models with melt transport, there is the porosity (which is transported with the melt velocity) and other compositional fields (like the depletion) that are transported with the solid velocity. At the moment, we can not prescribe different boundary conditions for different compositional fields. So we can only handle the outflow boundaries correctly for either the melt or the solid fields.

I don't really know what to do about this, except to enforce that 'Allow fixed composition on outflow boundaries' is set to true in models with melt and fix this later.

@tjhei
Copy link
Member

tjhei commented Jul 5, 2018

At the moment, we can not prescribe different boundary conditions for different compositional fields.

Yes, the problem is that we reuse the same sparsity pattern to avoid having to store n matrices for n compositional fields. One option would be to create the sparsity pattern only with the constraints that are relevant for all compositional fields and add the other constraints later on. This would waste a tiny amount of memory because most of the entries of those rows are zero, but I think we could care less in that situation...

@jdannberg
Copy link
Contributor Author

Okay, for now I disabled this in models with melt by default (because I don't have the time to implement what @tjhei suggested in the coming weeks, and fixing this in all other models is better than not fixing it at all).
It's not very pretty, so @gassmoeller let me know if there's a better solution.

@gassmoeller
Copy link
Member

@tjhei: for this PR your tester gives different results than the CIG tester, because master changed since this branch was created. We can only use the CIG tester results to update, because yours would not pass on master. Do you prefer it this way? We were confused by the differences for a few minutes. Do you want to change your tester configuration to test the merged version of the PR? The setting is under Branch Sources -> Github -> Discover pull requests from forks -> Strategy -> Merging the pull request with the current target branch revision.

@tjhei
Copy link
Member

tjhei commented Oct 3, 2018

This could be fixed be rebasing to the current master, right? This can be confusing, I admit. I set it up this way on purpose because it mirrors exactly what the old tester did before. I think there are pros and cons to both approaches, but let me think about it.

@tjhei
Copy link
Member

tjhei commented Oct 3, 2018

If we test pr-head then we might break master when merging a passing PR. But if we test pr-merge, things behave differently on your local machine and the tester, which can be confusing too. For example, the tester might not even compile while things work locally.

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a few minor comments left. Ready to merge after you addressed them.

void
Simulator<dim>::replace_outflow_boundary_ids(const unsigned int offset)
{
const QGauss<dim-1> quadrature_formula (finite_element.base_element(introspection.base_elements.temperature).degree+1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should be the velocity element, as you evaluate velocity below

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure. This is used to determine if the face has outflow when assembled for the temperature (or compositional) equation.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, disregard my comment then. There are identical by default anyway.

for (const auto row: current_constraints.get_local_lines())
{
if (current_constraints.is_constrained(row)
!= new_current_constraints.is_constrained(row))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think Juliane and Timo are right here. Apart from the case where no constraints change, we have to consider two cases:

  • The dof was constrained before and is not longer constrained in new_current_constraints. That means its relevant matrix entries will now be assembled and we have to ensure they are present. This is enforced by the current implementation.
  • The dof was not constrained before, and is now constrained in new_current_constraints. In this case the matrix entries exist, but they are not longer touched (as they exist in the constraint matrix), and therefore we can as well skip recreating the matrix (which is what is currently happening).

We did a test that worked although it should have triggered this bug (no_dirichlet_on_outflow_only_new_constraints) and so I am relatively confident the current implementation is correct.

current_constraints.close();
#endif

// let the melt handler add its constraints once before we solve the porosity system for the first time
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You moved or removed the following lines (because you do not need them any more?), please also remove the comment.

FEFaceValues<dim> fe_face_values (*mapping,
finite_element,
quadrature_formula,
update_gradients | update_values |
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you do not use the gradients, please remove the update flag

if (face->at_boundary())
{
Assert(face->boundary_id() <= offset,
ExcMessage("You can not use more than 10000 boundary indicators!"));
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think 10000 is not the correct number, somewhere above you write 128?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, the error should be something like "To allow 'Allow fixed composition on outflow boundaries' you are only allowed to use boundary ids between 0 and ".

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, the error should be "If you do not allow fixed composition on outflow boundaries you are only allowed...' I'll update it.


setup_system_matrix (introspection.index_sets.system_partitioning);
setup_system_preconditioner (introspection.index_sets.system_partitioning);
this->compute_current_constraints();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

you do not need this-> in this line.

@@ -0,0 +1,178 @@
# This test checks if fixed composition boundary conditions are only
# prescribed on outflow boundaries (if this option is activated in
# the input file), and in particular if this also works foe time-
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

foe -> for

@gassmoeller
Copy link
Member

@tjhei: You are right there are good reasons for both. I guess in the end we should just do both, and keep the difference in mind. When we find a way to circumvent the DNS issues with the CIG tester we can activate both test setups, which should make the differences more obvious (if they give different results).

@gassmoeller
Copy link
Member

Ok, looks like all comments have been addressed. Ready to merge.

@gassmoeller
Copy link
Member

@tjhei: Are you fine with the changes as well?

@tjhei tjhei merged commit b729283 into geodynamics:master Oct 6, 2018
freddrichards pushed a commit to freddrichards/aspect that referenced this pull request May 20, 2019
…tflow

[Bugfix] add an option that fixes temperature/composition only on inflow/close…
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants