From 265dfbea930a67ee952d8a5b38fbfa1d6c81bcea Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Tue, 21 Dec 2021 00:28:57 +0100 Subject: [PATCH] Allow inhomogeneous Dirichlet bcs,rhs and fix results.dox in step-60 --- doc/news/changes/minor/20211224MarcoFeder | 3 + examples/step-60/doc/intro.dox | 20 +-- examples/step-60/doc/results.dox | 169 ++++++++++++++++++---- examples/step-60/step-60.cc | 55 ++++--- 4 files changed, 188 insertions(+), 59 deletions(-) create mode 100644 doc/news/changes/minor/20211224MarcoFeder diff --git a/doc/news/changes/minor/20211224MarcoFeder b/doc/news/changes/minor/20211224MarcoFeder new file mode 100644 index 0000000000..6e0357de85 --- /dev/null +++ b/doc/news/changes/minor/20211224MarcoFeder @@ -0,0 +1,3 @@ +Changed: The step-60 tutorial now works also with non-zero Dirichlet BCs and non-zero right-hand side +
+(Marco Feder, 2021/12/24) diff --git a/examples/step-60/doc/intro.dox b/examples/step-60/doc/intro.dox index 45a2592e5c..23d3666bbe 100644 --- a/examples/step-60/doc/intro.dox +++ b/examples/step-60/doc/intro.dox @@ -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 boundary condition +$u_D$ on $\partial \Omega$, find the solution $u$ to @f{eqnarray*}{ -- \Delta u + \gamma^T \lambda &=& 0 \text{ in } \Omega\\ +- \Delta u + \gamma^T \lambda &=& f \text{ in } \Omega\\ \gamma u &=& g \text{ in } \Gamma \\ -u & = & 0 \text{ on } \partial\Omega. +u & = & u_D \text{ on } \partial\Omega. @f} -This is a constrained problem, where we are looking for a harmonic function $u$ -that satisfies homogeneous boundary conditions on $\partial\Omega$, subject to -the constraint $\gamma u = g$ using a Lagrange multiplier. +This is a constrained problem, where we are looking for a function $u$ that solves the +Poisson equation and that satisfies Dirichlet boundary conditions $u=u_D$ on $\partial \Omega$, +subject to the constraint $\gamma u = g$ using a Lagrange multiplier. -This problem has a physical interpretation: harmonic functions, i.e., functions +When $f=0$ this problem has a physical interpretation: harmonic functions, i.e., functions that satisfy the Laplace equation, can be thought of as the displacements of a membrane whose boundary values are prescribed. The current situation then corresponds to finding the shape of a membrane for which not only the @@ -98,7 +99,7 @@ conditions on $\partial\Omega$, we obtain the following variational problem: Given a sufficiently regular function $g$ on $\Gamma$, find the solution $u$ to @f{eqnarray*}{ -(\nabla u, \nabla v)_{\Omega} + (\lambda, \gamma v)_{\Gamma} &=& 0 \qquad \forall v \in V(\Omega) \\ +(\nabla u, \nabla v)_{\Omega} + (\lambda, \gamma v)_{\Gamma} &=& (f,v)_{\Omega} \qquad \forall v \in V(\Omega) \\ (\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} \qquad \forall q \in Q(\Gamma), @f} @@ -200,7 +201,7 @@ u \\ \end{pmatrix} = \begin{pmatrix} -0 \\ +F \\ G \end{pmatrix} \f] @@ -210,6 +211,7 @@ where @f{eqnarray*}{ K_{ij} &\dealcoloneq& (\nabla v_j, \nabla v_i)_\Omega \qquad i,j=1,\dots,n \\ C_{\alpha j} &\dealcoloneq& (v_j, q_\alpha)_\Gamma \qquad j=1,\dots,n, \alpha = 1,\dots, m \\\\ +F_{i} &\dealcoloneq& (f, v_i)_\Omega \qquad i=1,\dots,n \\ G_{\alpha} &\dealcoloneq& (g, q_\alpha)_\Gamma \qquad \alpha = 1,\dots, m. @f} diff --git a/examples/step-60/doc/results.dox b/examples/step-60/doc/results.dox index ae9c27c51c..1b4edef822 100644 --- a/examples/step-60/doc/results.dox +++ b/examples/step-60/doc/results.dox @@ -32,11 +32,11 @@ parameter file, we see the following: # --------------------- subsection Distributed Lagrange<1,2> set Coupling quadrature order = 3 + set Dirichlet boundary ids = 0, 1, 2, 3 set Embedded configuration finite element degree = 1 set Embedded space finite element degree = 1 set Embedding space finite element degree = 1 - set Homogeneous Dirichlet boundary ids = 0, 1, 2, 3 - set Initial embedded space refinement = 7 + set Initial embedded space refinement = 8 set Initial embedding space refinement = 4 set Local refinements steps near embedded domain = 3 set Use displacement in embedded interface = false @@ -125,6 +125,88 @@ subsection Distributed Lagrange<1,2> set Variable names = x,y,t end + subsection Embedding Dirichlet boundary conditions + # Sometimes it is convenient to use symbolic constants in the expression + # that describes the function, rather than having to use its numeric value + # everywhere the constant appears. These values can be defined using this + # parameter, in the form `var1=value1, var2=value2, ...'. + # + # A typical example would be to set this runtime parameter to + # `pi=3.1415926536' and then use `pi' in the expression of the actual + # formula. (That said, for convenience this class actually defines both + # `pi' and `Pi' by default, but you get the idea.) + set Function constants = + + # The formula that denotes the function you want to evaluate for + # particular values of the independent variables. This expression may + # contain any of the usual operations such as addition or multiplication, + # as well as all of the common functions such as `sin' or `cos'. In + # addition, it may contain expressions like `if(x>0, 1, -1)' where the + # expression evaluates to the second argument if the first argument is + # true, and to the third argument otherwise. For a full overview of + # possible expressions accepted see the documentation of the muparser + # library at http://muparser.beltoforion.de/. + # + # If the function you are describing represents a vector-valued function + # with multiple components, then separate the expressions for individual + # components by a semicolon. + set Function expression = 0 + + # The names of the variables as they will be used in the function, + # separated by commas. By default, the names of variables at which the + # function will be evaluated are `x' (in 1d), `x,y' (in 2d) or `x,y,z' (in + # 3d) for spatial coordinates and `t' for time. You can then use these + # variable names in your function expression and they will be replaced by + # the values of these variables at which the function is currently + # evaluated. However, you can also choose a different set of names for the + # independent variables at which to evaluate your function expression. For + # example, if you work in spherical coordinates, you may wish to set this + # input parameter to `r,phi,theta,t' and then use these variable names in + # your function expression. + set Variable names = x,y,t + end + + subsection Embedding rhs function + # Sometimes it is convenient to use symbolic constants in the expression + # that describes the function, rather than having to use its numeric value + # everywhere the constant appears. These values can be defined using this + # parameter, in the form `var1=value1, var2=value2, ...'. + # + # A typical example would be to set this runtime parameter to + # `pi=3.1415926536' and then use `pi' in the expression of the actual + # formula. (That said, for convenience this class actually defines both + # `pi' and `Pi' by default, but you get the idea.) + set Function constants = + + # The formula that denotes the function you want to evaluate for + # particular values of the independent variables. This expression may + # contain any of the usual operations such as addition or multiplication, + # as well as all of the common functions such as `sin' or `cos'. In + # addition, it may contain expressions like `if(x>0, 1, -1)' where the + # expression evaluates to the second argument if the first argument is + # true, and to the third argument otherwise. For a full overview of + # possible expressions accepted see the documentation of the muparser + # library at http://muparser.beltoforion.de/. + # + # If the function you are describing represents a vector-valued function + # with multiple components, then separate the expressions for individual + # components by a semicolon. + set Function expression = 0 + + # The names of the variables as they will be used in the function, + # separated by commas. By default, the names of variables at which the + # function will be evaluated are `x' (in 1d), `x,y' (in 2d) or `x,y,z' (in + # 3d) for spatial coordinates and `t' for time. You can then use these + # variable names in your function expression and they will be replaced by + # the values of these variables at which the function is currently + # evaluated. However, you can also choose a different set of names for the + # independent variables at which to evaluate your function expression. For + # example, if you work in spherical coordinates, you may wish to set this + # input parameter to `r,phi,theta,t' and then use these variable names in + # your function expression. + set Variable names = x,y,t + end + subsection Schur solver control set Log frequency = 1 set Log history = false @@ -141,15 +223,13 @@ If you now run the program, you will get a file called `used_parameters.prm`, containing a shorter version of the above parameters (without comments and documentation), documenting all parameters that were used to run your program: @code -# Parameter file generated with -# DEAL_II_PACKAGE_VERSION = 9.0.0 subsection Distributed Lagrange<1,2> set Coupling quadrature order = 3 + set Dirichlet boundary ids = 0, 1, 2, 3 set Embedded configuration finite element degree = 1 set Embedded space finite element degree = 1 set Embedding space finite element degree = 1 - set Homogeneous Dirichlet boundary ids = 0, 1, 2, 3 - set Initial embedded space refinement = 7 + set Initial embedded space refinement = 8 set Initial embedding space refinement = 4 set Local refinements steps near embedded domain = 3 set Use displacement in embedded interface = false @@ -164,6 +244,16 @@ subsection Distributed Lagrange<1,2> set Function expression = 1 set Variable names = x,y,t end + subsection Embedding Dirichlet boundary conditions + set Function constants = + set Function expression = 0 + set Variable names = x,y,t + end + subsection Embedding rhs function + set Function constants = + set Function expression = 0 + set Variable names = x,y,t + end subsection Schur solver control set Log frequency = 1 set Log history = false @@ -184,7 +274,7 @@ For example, you could use the following (perfectly valid) parameter file with this tutorial program: @code subsection Distributed Lagrange<1,2> - set Initial embedded space refinement = 7 + set Initial embedded space refinement = 8 set Initial embedding space refinement = 4 set Local refinements steps near embedded domain = 3 subsection Embedded configuration @@ -197,6 +287,16 @@ subsection Distributed Lagrange<1,2> set Function expression = 1 set Variable names = x,y,t end + subsection Embedding Dirichlet boundary conditions + set Function constants = + set Function expression = 0 + set Variable names = x,y,t + end + subsection Embedding rhs function + set Function constants = + set Function expression = 0 + set Variable names = x,y,t + end end @endcode @@ -230,55 +330,50 @@ as boundary of the portion of $\Omega$ inside $\Gamma$. Similarly on $\partial The output of the program will look like the following: @code -DEAL::Embedded dofs: 129 -DEAL::Embedding minimal diameter: 0.0110485, embedded maximal diameter: 0.00781250, ratio: 0.707107 +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.166266 -DEAL:cg::Convergence step 108 value 7.65958e-13 +DEAL:cg::Starting value 0.117692 +DEAL:cg::Convergence step 594 value 8.06558e-13 +---------------------------------------------+------------+------------+ -| Total CPU time elapsed since start | 0.586s | | +| Total CPU time elapsed since start | 1.48s | | | | | | | Section | no. calls | CPU time | % of total | +---------------------------------+-----------+------------+------------+ -| Assemble coupling system | 1 | 0.132s | 23% | -| Assemble system | 1 | 0.0733s | 12% | -| Output results | 1 | 0.087s | 15% | -| Setup coupling | 1 | 0.0244s | 4.2% | -| Setup grids and dofs | 1 | 0.0907s | 15% | -| Solve system | 1 | 0.178s | 30% | +| Assemble coupling system | 1 | 0.0381s | 2.6% | +| Assemble system | 1 | 0.153s | 10% | +| Output results | 1 | 0.11s | 7.4% | +| Setup coupling | 1 | 0.0364s | 2.5% | +| Setup grids and dofs | 1 | 0.168s | 11% | +| Solve system | 1 | 0.974s | 66% | +---------------------------------+-----------+------------+------------+ +---------------------------------------------+------------+------------+ -| Total wallclock time elapsed since start | 0.301s | | +| Total wallclock time elapsed since start | 0.798s | | | | | | | Section | no. calls | wall time | % of total | +---------------------------------+-----------+------------+------------+ -| Assemble coupling system | 1 | 0.0385s | 13% | -| Assemble system | 1 | 0.0131s | 4.3% | -| Output results | 1 | 0.0736s | 24% | -| Setup coupling | 1 | 0.0234s | 7.7% | -| Setup grids and dofs | 1 | 0.0679s | 23% | -| Solve system | 1 | 0.0832s | 28% | +| Assemble coupling system | 1 | 0.0469s | 5.9% | +| Assemble system | 1 | 0.0348s | 4.4% | +| Output results | 1 | 0.0821s | 10% | +| Setup coupling | 1 | 0.0371s | 4.7% | +| Setup grids and dofs | 1 | 0.157s | 20% | +| Solve system | 1 | 0.436s | 55% | +---------------------------------+-----------+------------+------------+ @endcode -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. - If the problem was set in a three-dimensional setting, and the immersed mesh was time dependent, it would be much more expensive to recreate the mesh at each step rather than use the technique we present here. Moreover, you may be able to create a very fast and optimized solver on a uniformly refined square or cubic grid, and embed the domain where you want to perform your computation using the technique presented here. This would require you to only have a surface -representatio of your domain (a much cheaper and easier mesh to produce). +representation of your domain (a much cheaper and easier mesh to produce). To play around a little bit, we are going to complicate a little the fictitious domain as well as the boundary conditions we impose on it. @@ -308,11 +403,21 @@ subsection Distributed Lagrange<1,2> set Function expression = x-.5 set Variable names = x,y,t end + subsection Embedding Dirichlet boundary conditions + set Function constants = + set Function expression = 0 + set Variable names = x,y,t + end + subsection Embedding rhs function + set Function constants = + set Function expression = 0.0 + set Variable names = x,y,t + end subsection Schur solver control set Log frequency = 1 set Log history = false set Log result = true - set Max steps = 100000 + set Max steps = 1000 set Reduction = 1.e-12 set Tolerance = 1.e-12 end diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc index 162f876f5c..7e28606448 100644 --- a/examples/step-60/step-60.cc +++ b/examples/step-60/step-60.cc @@ -247,12 +247,12 @@ namespace Step60 // $\Gamma$. unsigned int initial_embedded_refinement = 8; - // The list of boundary ids where we impose homogeneous Dirichlet boundary - // conditions. On the remaining boundary ids (if any), we impose - // homogeneous Neumann boundary conditions. - // As a default problem we have zero Dirichlet boundary conditions on + // The list of boundary ids where we impose (possibly inhomogeneous) + // Dirichlet boundary conditions. On the remaining boundary ids (if any), + // we impose homogeneous Neumann boundary conditions. As a default problem + // we have zero Dirichlet boundary conditions on // $\partial \Omega$ - std::list homogeneous_dirichlet_ids{0, 1, 2, 3}; + std::list dirichlet_ids{0, 1, 2, 3}; // FiniteElement degree of the embedding space: $V_h(\Omega)$ unsigned int embedding_space_finite_element_degree = 1; @@ -364,12 +364,22 @@ namespace Step60 std::unique_ptr> embedded_mapping; - // We do the same thing to specify the value of the function $g$, - // which is what we want our solution to be in the embedded space. + // We do the same thing to specify the value of the forcing term $f$. // In this case the Function is a scalar one. + ParameterAcceptorProxy> + embedding_rhs_function; + + // Then, we specify the value of the function $g$, which is what we want + // our solution to be in the embedded space. + // Again, in this case the Function is a scalar one. ParameterAcceptorProxy> embedded_value_function; + // Finally, the value of the Dirichlet boundary conditions on $\partial + // \Omega$ is specified. + ParameterAcceptorProxy> + embedding_dirichlet_boundary_function; + // Similarly to what we have done with the Functions::ParsedFunction class, // we repeat the same for the ReductionControl class, allowing us to // specify all possible stopping criteria for the Schur complement @@ -387,7 +397,7 @@ namespace Step60 AffineConstraints constraints; Vector solution; - Vector rhs; + Vector embedding_rhs; Vector lambda; Vector embedded_rhs; @@ -526,8 +536,7 @@ namespace Step60 add_parameter("Local refinements steps near embedded domain", delta_refinement); - add_parameter("Homogeneous Dirichlet boundary ids", - homogeneous_dirichlet_ids); + add_parameter("Dirichlet boundary ids", dirichlet_ids); add_parameter("Use displacement in embedded interface", use_displacement); @@ -556,7 +565,10 @@ namespace Step60 const Parameters ¶meters) : parameters(parameters) , embedded_configuration_function("Embedded configuration", spacedim) + , embedding_rhs_function("Embedding rhs function") , embedded_value_function("Embedded value") + , embedding_dirichlet_boundary_function( + "Embedding Dirichlet boundary conditions") , schur_solver_control("Schur solver control") , monitor(std::cout, TimerOutput::summary, TimerOutput::cpu_and_wall_times) { @@ -581,9 +593,15 @@ namespace Step60 "R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy"); }); + embedding_rhs_function.declare_parameters_call_back.connect( + []() -> void { ParameterAcceptor::prm.set("Function expression", "0"); }); + embedded_value_function.declare_parameters_call_back.connect( []() -> void { ParameterAcceptor::prm.set("Function expression", "1"); }); + embedding_dirichlet_boundary_function.declare_parameters_call_back.connect( + []() -> void { ParameterAcceptor::prm.set("Function expression", "0"); }); + schur_solver_control.declare_parameters_call_back.connect([]() -> void { ParameterAcceptor::prm.set("Max steps", "1000"); ParameterAcceptor::prm.set("Reduction", "1.e-12"); @@ -828,8 +846,7 @@ namespace Step60 // We now set up the DoFs of $\Omega$ and $\Gamma$: since they are // fundamentally independent (except for the fact that $\Omega$'s mesh is more - // refined "around" - // $\Gamma$) the procedure is standard. + // refined "around" $\Gamma$) the procedure is standard. template void DistributedLagrangeProblem::setup_embedding_dofs() { @@ -839,10 +856,10 @@ 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(), constraints); + *space_dh, id, embedding_dirichlet_boundary_function, constraints); } constraints.close(); @@ -852,7 +869,7 @@ namespace Step60 stiffness_sparsity.copy_from(dsp); stiffness_matrix.reinit(stiffness_sparsity); solution.reinit(space_dh->n_dofs()); - rhs.reinit(space_dh->n_dofs()); + embedding_rhs.reinit(space_dh->n_dofs()); deallog << "Embedding dofs: " << space_dh->n_dofs() << std::endl; } @@ -914,11 +931,13 @@ namespace Step60 { TimerOutput::Scope timer_section(monitor, "Assemble system"); - // Embedding stiffness matrix $K$, and the right hand side $G$. + // Embedding stiffness matrix $K$ and right hand sides $F$ and $G$. MatrixTools::create_laplace_matrix( *space_dh, QGauss(2 * space_fe->degree + 1), stiffness_matrix, + embedding_rhs_function, + embedding_rhs, static_cast *>(nullptr), constraints); @@ -978,9 +997,9 @@ namespace Step60 SolverCG> solver_cg(schur_solver_control); auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity()); - lambda = S_inv * embedded_rhs; + lambda = S_inv * (C * K_inv * embedding_rhs - embedded_rhs); - solution = K_inv * Ct * lambda; + solution = K_inv * (embedding_rhs - Ct * lambda); constraints.distribute(solution); } -- 2.39.5