From: Luca Heltai Date: Wed, 2 May 2018 23:09:29 +0000 (+0200) Subject: Added possibility to read parameter from command line. X-Git-Tag: v9.1.0-rc1~1190^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=188431486dd988ed0c01756dcc8e9eb5734b6c8c;p=dealii.git Added possibility to read parameter from command line. --- diff --git a/examples/step-60/doc/intro.dox b/examples/step-60/doc/intro.dox index 32443edac5..3fb137029e 100644 --- a/examples/step-60/doc/intro.dox +++ b/examples/step-60/doc/intro.dox @@ -41,7 +41,9 @@ The same is true, with a less regular range space (namely $H^{1/2}(\Gamma)$), when the dimension of $\Gamma$ is one less with respect to $\Omega$, and $\Gamma$ does not have a boundary. In this second case, the operator $\gamma$ is also known as the *trace* operator, and it is well defined for Lipschitz -co-dimension one curves and surfaces $\Gamma$ embedded in $\Omega$. +co-dimension one curves and surfaces $\Gamma$ embedded in $\Omega$ (read on +this wikipedia article +for further details on the trace operator). The co-dimension two case is a little more complicated, and in general it is not possible to construct a continous trace operator, not even from @@ -66,6 +68,19 @@ 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 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 +displacement at the boundary, but also on $\Gamma$ is prescribed. +For example, if $\Gamma$ is a closed curve in 2d space, then that would model +a soap film that is held in place by a wire loop along $\partial \Omega$ as +well as a second loop along $\Gamma$. In cases where $\Gamma$ is a whole area, +you can think of this as a membrane that is stretched over an obstacle where +$\Gamma$ is the contact area. (If the contact area is not known we have a +different problem -- called the "obstacle problem" -- which is modeled in +step-41.) + 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. @@ -84,7 +99,7 @@ Given a sufficiently regular function $g$ on $\Gamma$, find the solution $u$ to (\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} & \forall q \in Q(\Gamma), @f} -where $(\cdot, \cdot){\Omega}$ and $(\cdot, \cdot){\Gamma}$ represent, +where $(\cdot, \cdot)_{\Omega}$ and $(\cdot, \cdot)_{\Gamma}$ represent, respectively, $L^2$ scalar products in $\Omega$ and in $\Gamma$. Inspection of the variational formulation tells us that the space $V(\Omega)$ @@ -95,16 +110,23 @@ should be taken as $H^{1/2}(\Gamma)$. 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)$. +either $H^{-1}(\Gamma)$ or $H^{-1/2}(\Gamma)$. There are two options for the discretization of the problem above. One could choose matching discretizations, where the Triangulation for $\Gamma$ is aligned with the Triangulation for $\Omega$, or one could choose to discretize the two domains in a completely independent way. -While the first option is clearly more indicated for the simple problem we -proposed above, if the domain $\Gamma$ was to be time dependent, then the -second option could be a more viable solution. +The first option is clearly more indicated for the simple problem we +proposed above: it is sufficient to use a single Triangulation for $\Omega$ and +then impose certain constraints depending $\Gamma$. An example of this approach +is studied in step-40, where the solution has to stay above an obstacle and this +is achieved imposing constraints on $\Omega$. + +To solve more complex problems, for example one where the domain $\Gamma$ is time +dependent, the second option could be a more viable solution. Handling +non aligned meshes is complex by itself: to illustrate how is done we study a +simple problem. The technique we describe here is presented in the literature using one of many names: the **immersed finite element method**, the **fictitious boundary method**, the @@ -129,8 +151,8 @@ for any point $x \in \Gamma_0$, its coordinate $\psi(x)$ in $\Omega$; representing for any point $x$ the displacement vector applied in order to deform $x$ to its actual configuration $\psi(x) = x +\delta\psi(x)$. -We define the embedded reference domain $\Gamma_0$ `embedded_grid`, and on -this domain, we construct a finite dimensional space (`embedded_configuration_dh`) +We define the embedded reference domain $\Gamma_0$ `embedded_grid`: on +this triangulation we construct a finite dimensional space (`embedded_configuration_dh`) to describe either the deformation or the displacement through a FiniteElement system of FE_Q objects (`embedded_configuration_fe`). This finite dimensional space is used only to interpolate a user supplied function @@ -146,7 +168,7 @@ or a MappingQEulerian object are initialized with the `embedded_configuration` vector. In the embedding space, a standard finite dimensional space `space_dh` is -constructed on the embedding grid `space_grid` (`space_dh`), using the +constructed on the embedding grid `space_grid`, using the FiniteElement `space_fe`, following almost verbatim the approach taken in step-6. We represent the discretizations of the spaces $V$ and $Q$ with @@ -160,8 +182,9 @@ Q_h(\Gamma) = \text{span} \{q_i\}_{i=1}^m respectively, where $n$ is the dimension of `space_dh`, and $m$ the dimension of `embedded_dh`. -Once all the finite dimensional spaces are defined, we are left with the following -finite dimensional system of equations: +Once all the finite dimensional spaces are defined, the variational formulation +of the problem above leaves us with the following finite dimensional system +of equations: \f[ \begin{pmatrix} @@ -194,11 +217,11 @@ 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 mapping from $\hat K$ to $K$, and compute the integral -on $\hat K$ using a quadrature formula: +$C$, is computed on $\Gamma$. As usual in finite elements we split this +integral into contributions from all cells 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} @@ -213,7 +236,8 @@ 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): +To evaluate $(v_j \circ F_{K}) (\hat x_i)$ the following steps needs to be +taken (as shown in 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 @@ -226,21 +250,26 @@ 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: +The three steps above can be computed using 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 + cell + +- construct an FEValues object, with the given quadrature formula, and initialized + with the cell obtained in the first step. -- construct an FEValues object, with the given quadrature formula, and initialized with the cell obtained in the first step. +This is what the deal.II function VectorTools::point_value does when evaluating a +finite element field (not just a single shape function) at an arbitrary point; but +this is would be inefficient in this case. -The deal.II library offers a convenient wrapper to perform the first three steps +A better solution is to use 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 diff --git a/examples/step-60/doc/results.dox b/examples/step-60/doc/results.dox index ddaafdff0c..d582d0ac32 100644 --- a/examples/step-60/doc/results.dox +++ b/examples/step-60/doc/results.dox @@ -1,7 +1,8 @@

Results

-The first time the code is run, an exception is thrown and nothing is actually -computed. You'll get the following output: +Currently there is no `parameters.prm` file: objects can't be initialized. +For this reason 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: @@ -321,7 +322,7 @@ such that we have set Variable names = x,y,t end @endcode -produce the saddle on the right. +produces the saddle on the right. @@ -368,16 +369,23 @@ 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. +

Preconditioner

-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! +At the moment, we have no preconditioner on the Schur complement. This is ok for +two dimensional problems, where a few hundred iterations bring the residual down +to the machine precision, but it's not going to work in three dimensions. + +It is not obvious what a good preconditioner would be here. The physical problem +we are solving with the Schur complement, is to associate to the Dirichlet data +$g$, the value of the Lagrange multiplier $\lambda$. $\lambda$ can be +interpreted as the *jump* in the normal gradient that need to be imposed on $u$ +across $\Gamma$, in order to obtain the Dirichlet data $g$. + +So $S$ is some sort of Neumann to Dirichlet map, and we would like to have a +good approximation for the Dirichlet to Neumann map. A possibility would be to +use a Boundary Element approximation of the problem on $\Gamma$, and construct a +rough approximation of the hyper-singular operator for the Poisson problem +associated to $\Gamma$, which is precisely a Dirichlet to Neumann map.

Parallel Code

diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc index c3c24f016c..404e1444fc 100644 --- a/examples/step-60/step-60.cc +++ b/examples/step-60/step-60.cc @@ -398,11 +398,9 @@ namespace Step60 SparsityPattern stiffness_sparsity; SparsityPattern coupling_sparsity; - SparsityPattern embedded_sparsity; SparseMatrix stiffness_matrix; SparseMatrix coupling_matrix; - SparseMatrix embedded_stiffness_matrix; ConstraintMatrix constraints; @@ -876,10 +874,6 @@ namespace Step60 (parameters.embedded_space_finite_element_degree); embedded_dh->distribute_dofs(*embedded_fe); - DynamicSparsityPattern dsp(embedded_dh->n_dofs(), embedded_dh->n_dofs()); - DoFTools::make_sparsity_pattern(*embedded_dh, dsp); - embedded_sparsity.copy_from(dsp); - embedded_stiffness_matrix.reinit(embedded_sparsity); // By definition the rhs of the system we're solving involves only a zero // vector and $G$, which is computed using only $\Gamma$'s DoFs lambda.reinit(embedded_dh->n_dofs()); @@ -924,24 +918,23 @@ namespace Step60 { TimerOutput::Scope timer_section(monitor, "Assemble system"); - // Embedding stiffness matrix $K$ + // Embedding stiffness matrix $K$, and the right hand side $G$. MatrixTools::create_laplace_matrix(*space_dh, QGauss(2*space_fe->degree+1), stiffness_matrix, (const Function *) nullptr, constraints); - // Embedded stiffness matrix and rhs vector $G$ - MatrixTools::create_laplace_matrix(*embedded_mapping, - *embedded_dh, - QGauss(2*embedded_fe->degree+1), - embedded_stiffness_matrix, - embedded_value_function, - embedded_rhs); + VectorTools::create_right_hand_side(*embedded_mapping, + *embedded_dh, + QGauss(2*embedded_fe->degree+1), + embedded_value_function, + embedded_rhs); } { TimerOutput::Scope timer_section(monitor, "Assemble coupling system"); - // To compute the coupling matrix we use the NonMatching::create_coupling_mass_matrix - // tool, which works similarly to NonMatching::create_coupling_sparsity_pattern, - // requiring only an additional parameter: a constraint matrix + // To compute the coupling matrix we use the + // NonMatching::create_coupling_mass_matrix tool, which works similarly to + // NonMatching::create_coupling_sparsity_pattern, requiring only an + // additional parameter: a constraint matrix QGauss quad(parameters.coupling_quadrature_order); NonMatching::create_coupling_mass_matrix(*space_dh, *embedded_dh, @@ -968,7 +961,6 @@ namespace Step60 // Initializing the operators, as described in the introduction auto K = linear_operator(stiffness_matrix); - auto A = linear_operator(embedded_stiffness_matrix); auto Ct = linear_operator(coupling_matrix); auto C = transpose_operator(Ct); @@ -977,7 +969,7 @@ namespace Step60 // 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() ); + auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity()); lambda = S_inv * embedded_rhs; @@ -1023,7 +1015,7 @@ namespace Step60 -int main() +int main(int argc, char **argv) { try { @@ -1032,15 +1024,26 @@ int main() const unsigned int dim=1, spacedim=2; - // Differently to what happens in other tutorial programs, here we the + // Differently to what happens in other tutorial programs, here we use // ParameterAcceptor style of initialization, i.e., all objects are first // constructed, and then a single call to the static method // ParameterAcceptor::initialize is issued to fill all parameters of the // classes that are derived from ParameterAcceptor. + // + // We check if the user has specified a parameter file name to use when + // the program was launched. If so, try to read that parameter file, + // otherwise, try to read the file "parameters.prm". + // + // If the parameter file that was specified (implicitly or explicitly) + // does not exist, ParameterAcceptor::initialize will create one for you, + // and exit the program. DistributedLagrangeProblem::DistributedLagrangeProblemParameters parameters; DistributedLagrangeProblem problem(parameters); - ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm"); + std::string parameter_file = "parameters.prm"; + if (argc > 1) + parameter_file = argv[1]; + ParameterAcceptor::initialize(parameter_file, "used_parameters.prm"); problem.run(); } catch (std::exception &exc)