From 59192c0c41ac0b60ea2dd14e06f4f5fa85875d15 Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Wed, 2 May 2018 22:44:46 +0200 Subject: [PATCH] Written results, modified parameter file to use colorize --- examples/step-60/doc/intro.dox | 80 ++++-- examples/step-60/doc/results.dox | 409 +++++++++++++++++++++++++++++-- examples/step-60/step-60.cc | 26 +- 3 files changed, 457 insertions(+), 58 deletions(-) diff --git a/examples/step-60/doc/intro.dox b/examples/step-60/doc/intro.dox index 221f2b5c5a..32443edac5 100644 --- a/examples/step-60/doc/intro.dox +++ b/examples/step-60/doc/intro.dox @@ -66,6 +66,10 @@ 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. +As a first example we study the zero Dirichlet boundary condition on +$\partial\Omega$. The same equations apply if we apply zero Neumann boundary +conditions on $\partial\Omega$ or a mix of the two. + The variational formulation can be derived by introducing two infinite dimensional spaces $V(\Omega)$ and $Q^*(\Gamma)$, respectively for the solution $u$ and for the Lagrange multiplier $\lambda$. @@ -88,7 +92,7 @@ can be taken to be $H^1_0(\Omega)$. The space $Q(\Gamma)$, in the co-dimension zero case, should be taken as $H^1(\Gamma)$, while in the co-dimension one case should be taken as $H^{1/2}(\Gamma)$. -The function $g$ should, therefore, should be either in $H^1(\Gamma)$ (for the +The function $g$ should therefore be either in $H^1(\Gamma)$ (for the co-dimension zero case) or $H^{1/2}(\Gamma)$ (for the co-dimension one case). This leaves us with a Lagrange multiplier $\lambda$ in $Q^*(\Gamma)$, which is either $(H^1(\Gamma))^*$ or $H^{-1/2}(\Gamma)$. @@ -183,17 +187,18 @@ C_{\alpha j} &:=& (v_j, q_\alpha)_\Gamma &j=1,\dots,n, \alpha = 1,\dots, m \\\\ G_{\alpha} &:=& (g, q_\alpha)_\Gamma & \alpha = 1,\dots, m. @f} -While the matrix $K$ is the standard stiffness matrix for the Poisson problem -on $\Omega$, and the vector $G$ is a standard right-hand-side vector for a -finite element problem with forcing term $g$ on $\Gamma$, (see, for example, -step-3), the matrix $C$ or its transpose $C^T$ are non-standard since they -couple information on two non-matching grids. +While the matrix $K$ is the standard stiffness matrix for the Poisson problem on +$\Omega$, and the vector $G$ is a standard right-hand-side vector for a finite +element problem with forcing term $g$ on $\Gamma$, (see, for example, step-3), +the matrix $C$ or its transpose $C^T$ are non-standard since they couple +information on two non-matching grids. -In particular, the integral that appears in the computation of a single entry of $C$, -is computed on $\Gamma$. As usual in finite elements, we split this integral on each -cell of the triangulation used to discretize $\Gamma$, we tranform the integral on $K$ to -an integral on the reference element $\hat K$, where $F_{K}$ is the corresponding -shape function, and compute the integral there using a quadrature formula: +In particular, the integral that appears in the computation of a single entry of +$C$, is computed on $\Gamma$. As usual in finite elements, we split this +integral on each cell of the triangulation used to discretize $\Gamma$, we +tranform the integral on $K$ to an integral on the reference element $\hat K$, +where $F_{K}$ is the mapping from $\hat K$ to $K$, and compute the integral +on $\hat K$ using a quadrature formula: \f[ C_{\alpha j} := (v_j, q_\alpha)_\Gamma = \sum_{K\in \Gamma} \int_{\hat K} @@ -201,19 +206,48 @@ C_{\alpha j} := (v_j, q_\alpha)_\Gamma = \sum_{K\in \Gamma} \int_{\hat K} \sum_{K\in \Gamma} \sum_{i=1}^{n_q} \big(\hat q_\alpha(\hat x_i) (v_j \circ F_{K}) (\hat x_i) J_K (\hat x_i) w_i \big) \f] -Computing this sum is difficult because we have to evaluate $(v_j \circ F_{K}) (\hat x_i)$; this -complicated process is illustrated in the following picture: +Computing this sum is non-trivial because we have to evaluate $(v_j \circ F_{K}) +(\hat x_i)$. In general, if $\Gamma$ and $\Omega$ are not aligned, the point +$F_{K}(\hat x_i)$ is completely arbitrary with respect to $\Omega$, and unless +we figure out a way to interpolate all basis functions of $V_h(\Omega)$ on an +arbitrary point on $\Omega$, we cannot compute the integral needed for an entry +of the matrix $C$. + +The process is described by the following algorithm (and by the picture below): -

- -

+- For a given cell $K$ in $\Gamma$ compute the real point $y_i \coloneqq F_{K} +(\hat x_i)$, where $x_i$ is one of the quadrature points used for the integral +on $K \subseteq \Gamma$. -This translates in the following algorithm: -- Fix the cell $K$ in $\Gamma$ and compute the real point $y_i \coloneqq F_{K} (\hat x_i)$ - Find the cell of $\Omega$ in which $y_i$ lies, say it's $T$. -- To evaluate the basis function use invert shape function $G_{T}$ on $T$: - $v_j(y_i) = \hat v_j \circ G^{-1}_{T} (y_i)$. -We solve the above saddle point problem by iterating over the Schur complement -(which is described, for example, in step-20), -and we construct such complement using LinearOperator classes. +- To evaluate the basis function use the inverse of the mapping $G_T$ that +transforms the reference element $\hat T$ into the element $T$: $v_j(y_i) = \hat +v_j \circ G^{-1}_{T} (y_i)$. + +

+ +The three steps above can be computed by calling, in turn: + +- GridTools::find_active_cell_around_point, followed by + +- Mapping::tranform_real_to_unit_cell + +- construct a custom Quadrature formula, containing the point in the reference +cell + +- construct an FEValues object, with the given quadrature formula, and initialized with the cell obtained in the first step. + +The deal.II library offers a convenient wrapper to perform the first three steps +on a collection of points: GridTools::compute_point_locations. If one is +actually interested in computing the full coupling matrix, then it is possible +to call the method NonMatching::create_coupling_mass_matrix, that performs the +above steps in an efficient way, reusing all possible data structures, and +gathering expensive steps together. This is the function we'll be using later +in this tutorial. + +We solve the final saddle point problem by an iterative solver, applied to the +Schur complement $S$ (whose construction is described, for example, in step-20), and +we construct $S$ using LinearOperator classes. diff --git a/examples/step-60/doc/results.dox b/examples/step-60/doc/results.dox index c48f2a2d22..ddaafdff0c 100644 --- a/examples/step-60/doc/results.dox +++ b/examples/step-60/doc/results.dox @@ -1,26 +1,383 @@

Results

-

Test case 1:

- -For the default problem the value of u on Gamma is 1: this is like imposing a constant -Dirichlet boundary on the portion of $\Omega$ inside gamma. -

- -

-In this second image we can appreciate how, using GridTools::compute_point_locations , -we have been able to adaptively refine $\Omega$ exactly where the solution is varying -the most: - -

- -

+The first time the code is run, an exception is thrown and nothing is actually +computed. You'll get the following output: +@code +---------------------------------------------------- +Exception on processing: + +-------------------------------------------------------- +An error occurred in line <74> of file <../source/base/parameter_acceptor.cc> in function + static void dealii::ParameterAcceptor::initialize(const std::string &, const std::string &, const ParameterHandler::OutputStyle, dealii::ParameterHandler &) +The violated condition was: + false +Additional information: + You specified as input parameter file, but it does not exist. We created it for you. +-------------------------------------------------------- + +Aborting! +---------------------------------------------------- +@endcode + +By inspecting the parameter file, we see the following: + +@code +# Listing of Parameters +# --------------------- +subsection Distributed Lagrange<1,2> + set Coupling quadrature order = 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 embedding space refinement = 4 + set Local refinements steps near embedded domain = 3 + set Use displacement in embedded interface = false + set Verbosity level = 10 + + + subsection Embedded configuration + # 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 = R=.3, Cx=.4, Cy=.4 # default: + + # 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 = R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy # default: 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 Embedded value + # 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 = 1 # default: 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 + set Log result = true + set Max steps = 1000 # default: 100 + set Reduction = 1.e-12 # default: 1.e-2 + set Tolerance = 1.e-12 # default: 1.e-10 + end + +end +@endcode + +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-pre +subsection Distributed Lagrange<1,2> + set Coupling quadrature order = 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 embedding space refinement = 4 + set Local refinements steps near embedded domain = 3 + set Use displacement in embedded interface = false + set Verbosity level = 10 + subsection Embedded configuration + set Function constants = R=.3, Cx=.4, Cy=.4 + set Function expression = R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy + set Variable names = x,y,t + end + subsection Embedded value + set Function constants = + set Function expression = 1 + 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 = 1000 + set Reduction = 1.e-12 + set Tolerance = 1.e-12 + end +end +@endcode + +The rationale behind creating first `parameters.prm` file (the first time the +program is run) and then a `used_parameters.prm` (every other times you run the +program), is because you may want to leave most parameters to their default +values, and only modify a handful of them. + +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 embedding space refinement = 4 + set Local refinements steps near embedded domain = 3 + subsection Embedded configuration + set Function constants = R=.3, Cx=.4, Cy=.4 + set Function expression = R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy + set Variable names = x,y,t + end + subsection Embedded value + set Function constants = + set Function expression = 1 + set Variable names = x,y,t + end +end +@endcode + +and you would obtain exactly the same results as in test case 1 below. + +

Test case 1:

+ +For the default problem the value of $u$ on Gamma is set to the constant $1$: +this is like imposing a constant Dirichlet boundary condition on $\Gamma$, seen +as boundary of the portion of $\Omega$ inside $\Gamma$. Similarly on $\partial +\Omega$ we have zero Dirichlet boundary conditions. + + + + + + +
+ + + +
+ +The output of the program will look like the following on your pc: + +@code +DEAL::Embedded dofs: 129 +DEAL::Embedding minimal diameter: 0.0110485, embedded maximal diameter: 0.00781250, ratio: 0.707107 +DEAL::Embedding dofs: 2429 +DEAL:cg::Starting value 0.166266 +DEAL:cg::Convergence step 108 value 7.65958e-13 + + ++---------------------------------------------+------------+------------+ +| Total CPU time elapsed since start | 0.586s | | +| | | | +| 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% | ++---------------------------------+-----------+------------+------------+ + + + ++---------------------------------------------+------------+------------+ +| Total wallclock time elapsed since start | 0.301s | | +| | | | +| 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% | ++---------------------------------+-----------+------------+------------+ + +@endcode + +You may notice that, in terms of CPU time, assembling the coupling system is +twice more expensive than assembling the standard Poisson system, even though +the matrix is smaller. This is due to the non-matching nature of the +discretization. Wether 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). + +To play around a little bit, we are going to complicate a little the ficticious +domain as well as the boundary conditions we impose on it. + +

Test case 2 and 3:

+ +If we use the following parameter file: +@code +subsection Distributed Lagrange<1,2> + set Coupling quadrature order = 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 = 8 + set Initial embedding space refinement = 4 + set Local refinements steps near embedded domain = 4 + set Use displacement in embedded interface = false + set Verbosity level = 10 + subsection Embedded configuration + set Function constants = R=.3, Cx=.5, Cy=.5, r=.1, w=12 + set Function expression = (R+r*cos(w*pi*x))*cos(2*pi*x)+Cx; (R+r*cos(w*pi*x))*sin(2*pi*x)+Cy + set Variable names = x,y,t + end + subsection Embedded value + set Function constants = + set Function expression = x-.5 + 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 Reduction = 1.e-12 + set Tolerance = 1.e-12 + end +end +@endcode + +we get a "flowery" looking domain, where we impose a linear boundary condition +$g=x-.5$. This test shows that the method is actually quite accurate in +recovering an exactly linear function from its boundary conditions, and even +though the meshes are not aligned, we obtain a pretty good result. + +Replacing $x-.5$ with $2(x-.5)^2-2(y-.5)^2$, i.e., modifying the parameter file +such that we have +@code + ... + subsection Embedded value + set Function constants = + set Function expression = 2*(x-.5)^2-2*(y-.5)^2 + set Variable names = x,y,t + end +@endcode +produce the saddle on the right. + + + + + + +
+ + + +

Possibilities for extensions

-Add something +The `parameters.prm` file is a perfect tool to play with: we can easily impose a +Neumann boundary, change the parametrization, the function values and the +refinements. + +However, non aligned meshes are complex to handle, as you shall see after some +changes. For example set the following parameters: + +@code +subsection Distributed Lagrange<1,2> + set Homogeneous Dirichlet boundary ids = 0, 1, 3 + subsection Embedded configuration + set Function constants = R=.2, Cx=.5, Cy=.4 + set Function expression = R*sin(2*pi*x)+0.5*R*sin(4*pi*x)+Cx; 1.5*R*cos(2*pi*x)+Cy + end +end +@endcode + +Keeping the default refinement options you shall notice some glitch in the upper part +of this hot air balloon shape. Why? What's the difference with the following? + + + + + + +
+ + + +
+ +To solve this problem we need to improve the way in which adaptive refinement is +performed. Currently the vector of support points contains the vertices of +$\Gamma$, if the finite dimensional space for the configuration is of degree +one: these are too far away to make $\Omega$ refine properly and simply +increasing the refinements of $\Gamma$ will not be enough, the conjugate +gradient won't converge. + +As you modify the code you shall notice that, for instance, with complex +geometries there shall be other convergence problems: parameters are handy, but do +not solve all problems!

Parallel Code

@@ -33,14 +390,18 @@ parallel::distributed::Triangulation ). When using non-matching grids in parallel a problem arises: to compute the matrix $C$ a process needs information about both meshes on the same portion of real space but, when working with -distributed meshes, this might not be the case. +distributed meshes, this information may not be available. + Various strategies can be implemented to tackle this problem: - distribute the two meshes so that this constraint is satisfied -- use communication for the parts of real space where the constraint is - not satisfied -- make use of a shared triangulation and a distributed triangulation -The latter strategy is clearly the easier to implement, as all -the function used in this tutorial program can work letting $\Omega$ -be distributed and $\Gamma$ be a shared triangulation. +- use communication for the parts of real space where the constraint is not + satisfied + +- use a shared triangulation for the embedding space, and a distributed + triangulation for the emdedded configuration + +The latter strategy is clearly the easiest to implement, as most of the +functions used in this tutorial program will work unchanged also in the parallel +case. diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc index 2c05c8c39c..c3c24f016c 100644 --- a/examples/step-60/step-60.cc +++ b/examples/step-60/step-60.cc @@ -13,7 +13,7 @@ * * --------------------------------------------------------------------- * - * Author: Luca Heltai, Giovanni Alzetta, + * Authors: Luca Heltai, Giovanni Alzetta, * International School for Advanced Studies, Trieste, 2018 */ @@ -243,10 +243,12 @@ namespace Step60 // $\Gamma$. unsigned int initial_embedded_refinement = 7; - // A list of boundary ids where we impose homogeneous Dirichlet boundary + // 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 - std::list homogeneous_dirichlet_ids {0}; + // 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}; // FiniteElement degree of the embedding space: $V_h(\Omega)$ unsigned int embedding_space_finite_element_degree = 1; @@ -648,7 +650,10 @@ namespace Step60 // Initializing $\Omega$: // constructing the Triangulation and wrapping it into a unique_ptr space_grid = std_cxx14::make_unique >(); - GridGenerator::hyper_cube(*space_grid); + // The last argument is set to true: this activates colorize, which + // we use to assign the Dirichlet and Neumann conditions. See + // GridGenerator::hyper_rectangle for the details on the ids used. + GridGenerator::hyper_cube(*space_grid,0,1,true); // Requesting the varius values to the parameters object, which is // of type DistributedLagrangeProblemParameters space_grid->refine_global(parameters.initial_refinement); @@ -739,6 +744,9 @@ namespace Step60 *embedded_dh, support_points); + for (unsigned int i=1; in_dofs(); ++i) + support_points.emplace_back((support_points[i-1]+support_points[i])/2); + // Once we have the support points of the embedded finite element space, we // would like to identify what cells of the embedding space contain what // support point, to get a chance at refining the embedding grid where it is @@ -958,9 +966,6 @@ namespace Step60 SparseDirectUMFPACK K_inv_umfpack; K_inv_umfpack.initialize(stiffness_matrix); - // Same thing, for the embedded space -// SparseDirectUMFPACK A_inv_umfpack; -// A_inv_umfpack.initialize(embedded_stiffness_matrix); // Initializing the operators, as described in the introduction auto K = linear_operator(stiffness_matrix); auto A = linear_operator(embedded_stiffness_matrix); @@ -968,12 +973,11 @@ namespace Step60 auto C = transpose_operator(Ct); auto K_inv = linear_operator(K, K_inv_umfpack); -// auto A_inv = linear_operator(A, A_inv_umfpack); - auto S = C*K_inv*Ct; // Using the Schur complement method + auto S = C*K_inv*Ct; SolverCG > solver_cg(schur_solver_control); - auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity() );//A_inv); + auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity() ); lambda = S_inv * embedded_rhs; -- 2.39.5