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

Allow inhomogeneous Dirichlet in Step-60 #13105

Merged
merged 1 commit into from Dec 28, 2021

Conversation

fdrmrc
Copy link
Contributor

@fdrmrc fdrmrc commented Dec 21, 2021

This PR fixes #13034 and allows inhomogeneous Dirichlet bcs. Indeed, with the current version, the rhs vector is not used and changing only the call to VectorTools::interpolate_boundary_values() doesn't produce the right result as the inhomogeneous constraints are not distributed to the rhs. Now it works even with non-constant bcs, like the one in the picture.
sin_cos

Copy link
Member

@luca-heltai luca-heltai left a comment

Choose a reason for hiding this comment

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

I think you are still missing a small change in the assemble_system function, since you also need to pass the embedding_rhs to the MatrixTools::create_laplace_matrix function

@@ -125,6 +125,47 @@ subsection Distributed Lagrange<1,2>
set Variable names = x,y,t
end

subsection Embedding value
Copy link
Member

Choose a reason for hiding this comment

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

Maybe I would rename this to Embedding Dirichlet boundary conditions

@@ -308,6 +357,11 @@ subsection Distributed Lagrange<1,2>
set Function expression = x-.5
set Variable names = x,y,t
end
subsection Embedding value
Copy link
Member

Choose a reason for hiding this comment

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

Same here: Embedding Dirichlet boundary conditions

// In this case the Function is a scalar one.
ParameterAcceptorProxy<Functions::ParsedFunction<spacedim>>
embedding_value_function;
Copy link
Member

Choose a reason for hiding this comment

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

embedding_dirichlet_boundary_conditions

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That's related to the value of the forcing term f, the first component of the right-hand side, which is named embedding_rhs, not to the bcs.

Copy link
Member

Choose a reason for hiding this comment

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

If I understand this correctly, this is the rhs of the Poisson problem, right? I would not name this in the same way as the other term, simply because they have a different meaning. Maybe I would rename this embedding_rhs_function, and similarly for the parameter files.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Right, it's the rhs for Poisson. I'll follow your suggestion.

Comment on lines 589 to 598
embedding_value_function.declare_parameters_call_back.connect(
[]() -> void { ParameterAcceptor::prm.set("Function expression", "0"); });

Copy link
Member

Choose a reason for hiding this comment

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

This is not necessary, since by default the function is set to "0"

@@ -839,7 +850,7 @@ namespace Step60
space_dh->distribute_dofs(*space_fe);

DoFTools::make_hanging_node_constraints(*space_dh, constraints);
for (auto id : parameters.homogeneous_dirichlet_ids)
for (auto id : parameters.dirichlet_ids)
{
VectorTools::interpolate_boundary_values(
*space_dh, id, Functions::ZeroFunction<spacedim>(), constraints);
Copy link
Member

Choose a reason for hiding this comment

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

Shouldn't this be the embedding_dirichlet_boundary_values instead of Functions::ZeroFunction<spacedim>()?

Copy link
Contributor Author

@fdrmrc fdrmrc Dec 21, 2021

Choose a reason for hiding this comment

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

Here I was assuming that the user should create its own function object to specify different bcs. Anyway, you're totally right, I should create another parsed function for that. I'm fixing this

Copy link
Member

Choose a reason for hiding this comment

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

Ah, ok. In the comment to the PR you talked about Dirichlet BCs, but your changes actually add a non homogeneous rhs, rather than non-homogeneous Dirichlet BCs.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

If the user doesn't change Functions::ZeroFunction<spacedim>(), then you're right: I only allowed a possibly non-zero rhs.

But with my changes, if the "f" in the original problem is 0, and the user specifies a non-homogeneous Dirichlet BCs on dOmega, the call to create_laplace_matrix() will distribute correctly the inhomogeneous constraints to the rhs vector. Before this PR, those constraints were not distributed in this scenario. Of course everything works with a non-zero "f" too, like in the picture above.

@fdrmrc
Copy link
Contributor Author

fdrmrc commented Dec 21, 2021

Thanks Luca for the review.

I think you are still missing a small change in the assemble_system function, since you also need to pass the embedding_rhs to the MatrixTools::create_laplace_matrix function

Here, embedding_value_function is related to the first component of the right-hand side, namely embedding_rhs, and both are used in the call to MatrixTools::create_laplace_matrix().

@luca-heltai
Copy link
Member

Ciao Marco, welcome! Please, make sure you run "make indent" before committing. The continuous integration system is complaining that step-60 is not indented correctly. Other than this, this looks fine to me!

@fdrmrc
Copy link
Contributor Author

fdrmrc commented Dec 21, 2021

I noticed one last thing here:

You may notice that, in terms of CPU time, assembling the coupling system is
twice as expensive as assembling the standard Poisson system, even though the
matrix is smaller. This is due to the non-matching nature of the discretization.
Whether this is acceptable or not, depends on the applications.

This is not what I observe on my machine after several tests. For instance, the default .prm file gives:

DEAL::Embedded dofs: 257
DEAL::Embedding minimal diameter: 0.0110485, embedded maximal diameter: 0.00736292, ratio: 0.666416
DEAL::Embedding dofs: 2429
DEAL:cg::Starting value 0.117692
DEAL:cg::Convergence step 598 value 7.68984e-13


+---------------------------------------------+------------+------------+
| Total CPU time elapsed since start          |     0.989s |            |
|                                             |            |            |
| Section                         | no. calls |  CPU time  | % of total |
+---------------------------------+-----------+------------+------------+
| Assemble coupling system        |         1 |     0.037s |       3.7% |
| Assemble system                 |         1 |    0.0877s |       8.9% |
| Output results                  |         1 |     0.116s |        12% |
| Setup coupling                  |         1 |     0.032s |       3.2% |
| Setup grids and dofs            |         1 |     0.161s |        16% |
| Solve system                    |         1 |     0.551s |        56% |
+---------------------------------+-----------+------------+------------+

Can I delete these lines?

@fdrmrc
Copy link
Contributor Author

fdrmrc commented Dec 22, 2021

I've added the function for Dirichlet BCs. The only thing that has to be changed is the comment above about CPU times.

The indentation failure was coming from some whitespaces in results.dox file.

@luca-heltai
Copy link
Member

Could you address my last comment, and also change the documentation of the tutorial, including the description of the equations? We now have a f in the rhs, and a u_D as boundary data. Otherwise the program does not do what is described by the equations at the beginning.

@fdrmrc
Copy link
Contributor Author

fdrmrc commented Dec 23, 2021

Sure, thanks. Should I delete those lines about CPU times?

@luca-heltai
Copy link
Member

Yes, indeed, please adjust the comment about the CPU time. A lot of things have changed in the library since the implementation of step-60. :)

@fdrmrc
Copy link
Contributor Author

fdrmrc commented Dec 24, 2021

Everything should be okay now!

@luca-heltai
Copy link
Member

One last thing: can you squash all your commits into one, and add a file in the directory: doc/news/changes/minor with a short description of your change? You already have one file there, so you should now the drill. :)

I'll let someone else comment on this before merging, since you are one of my PhD students. :)

Copy link
Member

@jppelteret jppelteret left a comment

Choose a reason for hiding this comment

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

Looks good to me, thanks @fdrmrc . I have just a couple of very small questions, which shouldn't take long to resolve.

@@ -0,0 +1,3 @@
Changed: The step-60 tutorial now works also with non-zero Dirichlet BCs and non-zero right-hand side
<br>
(Marco Feder, 2021/12/24)
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 that you need a new line at the end of this file.

@@ -60,19 +60,20 @@ embedded domain (`dim`) is always smaller by one or equal with respect to the
dimension of the embedding domain $\Omega$ (`spacedim`).

We are going to solve the following differential problem: given a sufficiently
regular function $g$ on $\Gamma$, find the solution $u$ to
regular function $g$ on $\Gamma$, a forcing term $f \in L^2(\Omega)$ and a Dirichlet data
Copy link
Member

Choose a reason for hiding this comment

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

a Dirichlet data

This doesn't read right to me. Do you mean "Dirichlet conditions", or something else?

Copy link
Contributor Author

@fdrmrc fdrmrc Dec 27, 2021

Choose a reason for hiding this comment

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

Yes, exactly. I could replace that with "a Dirichlet boundary condition", since that is imposed on dOmega only. Does that work for you?

Copy link
Member

Choose a reason for hiding this comment

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

Yes, that would be perfect thanks!

@@ -125,6 +125,88 @@ subsection Distributed Lagrange<1,2>
set Variable names = x,y,t
end

subsection Embedding Dirichlet boundary conditions
Copy link
Member

Choose a reason for hiding this comment

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

The documentation/commentary for these new set functions is exactly the same as what's found elsewhere in the tutorials, right?

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

Copy link
Member

Choose a reason for hiding this comment

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

Great, thanks for the confirmation. I did a quick check and it looked identical, but I thought I'd just make sure.

Copy link
Member

Choose a reason for hiding this comment

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

Ah, I see that you did this already. Awesome :-)

Copy link
Member

@masterleinad masterleinad left a comment

Choose a reason for hiding this comment

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

Apart from @jppelteret comments, this looks good to me.

@fdrmrc fdrmrc force-pushed the Step-60_inhomogeneous_Dirichlet branch from 56d5aab to 265dfbe Compare December 28, 2021 15:15
@fdrmrc
Copy link
Contributor Author

fdrmrc commented Dec 28, 2021

Thanks for the checks, @masterleinad and @jppelteret. Everything should be okay now

Copy link
Member

@jppelteret jppelteret left a comment

Choose a reason for hiding this comment

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

Thanks a lot for the nice extension @fdrmrc!

@jppelteret
Copy link
Member

/rebuild

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Wrong .prm files in step-60
4 participants