Documentation of the Config Files for Shape Optimization Problems#
Let us take a detailed look at the config files for shape optimization problems and
discusss the corresponding parameters. The corresponding
config file used for this discussion is config.ini
,
which is the config file used for Shape Optimization with a Poisson Problem.
For shape optimization problems, the config file is a lot larger compared to the config files for optimal control. However, several important parameters are shared between both _typing of optimization problems.
A general config file for shape optimization has the following sections Mesh, StateSystem, OptimizationRoutine, AlgoLBFGS, AlgoCG, ShapeGradient, Regularization, MeshQuality and Output. We go over these sections and each parameter in them in the following.
As in Documentation of the Config Files for Optimal Control Problems, we refer to the documentation of the configparser module for a detailed description of how these config files can be structured. Moreover, we remark that cashocs has a default behavior for almost all of these parameters, which is triggered when they are NOT specified in the config file, and we will discuss this behavior for each parameter in this tutorial. For a summary over all parameters and their default values look at the end of this page.
Section Mesh#
The mesh section is, in contrast to the corresponding one for optimal control problems, more important and also has more parameters. However, its primary importance is for remeshing, which we cover in Remeshing with cashocs.
As first parameter, we have
mesh_file = ./mesh/mesh.xdmf
This specifies a path to a .xdmf file containing the discretized geometry. For all purposes, cashocs assumes that this .xdmf file was generated via conversion from a
GMSH file using the command line command cashocs.convert()
.
Note, that the corresponding files for the boundaries and subdomains are generated
automatically with cashocs.convert()
, and they will also be read by import_mesh
if they are present.
The second parameter in the Mesh section, gmsh_file
, is defined via
gmsh_file = ./mesh/mesh.msh
This defines the path to the GMSH .msh file which was used to create the .xdmf file
specified in mesh_file
. As before, this parameter is only relevant for remeshing
purposes, and not needed otherwise.
The next parameter is geo_file
, which is the final file we need for remeshing (
and only there). It is also given by a path to a file, in this case to the GMSH .geo
file used to generate the gmsh_file
. It is specified, .e.g., as
geo_file = ./mesh/mesh.geo
Note
For a detailed discussion of how to use these parameters we refer to Remeshing with cashocs.
Next up is a boolean flag that is used to indicate whether remeshing shall be performed
remesh = False
As the remeshing feature is experimental, we do advise to always try without
remeshing. Note, that by default this flag is set to remehs = False
so that remeshing is disabled.
Finally, we have the boolean flag show_gmsh_output
, specified via
show_gmsh_output = False
This is used to toggle on / off the terminal output of GMSH when it performs a
remeshing operation. This can be helpful for debugging purposes. By default, this
is set to show_gmsh_output = False
.
As stated throughout the Mesh section, these parameters are optional most of the time, and are only really required for remeshing. You can safely leave them out of your config file, and you should not need them, unless you want to perform remeshing.
Section StateSystem#
The StateSystem section is in complete analogy to the corresponding one for optimal control problems. For the sake of completeness, we briefly recall the parameters here, anyway.
The first parameter is is_linear
, and can be set as
is_linear = True
This is a boolean flag that indicates whether the state system is linear or not.
The default value for this parameter is is_linear = False
, as every linear problem can also be
interpreted as a nonlinear one.
The next parameters are used to define the tolerances of the Newton solver, in case a nonlinear state system has to be solved
newton_rtol = 1e11
newton_atol = 1e13
Here, newton_rtol
sets the relative, and newton_atol
the absolute tolerance
for Newton’s method. Their default values are newton_rtol = 1e11
and
newton_atol = 1e13
.
The next parameter for the Newton iteration is the maximum number of iterations it is allowed to perform before the iteration is cancelled. This is controlled via
newton_iter = 50
which defaults to newton_iter = 50
.
The parameter newton_damped
, which is set via
newton_damped = True
is a boolean flag, indicating whether a damping strategy should be performed for the
Newton method, or whether the classical NewtonRaphson iteration shall be used. This
defaults to newton_damped = False
(as this is faster), but for some problems it might be beneficial to
use damping in order to enhance the convergence of the nonlinear solver.
Additionally, we have the boolean parameter newton_inexact
, defined via
newton_inexact = False
which sets up an inexact Newton method for solving nonlinear problems in case this is newton_inexact = True
. The default is newton_inexact = False
.
Next, we have the parameter
newton_verbose = False
This is used to toggle the verbose output of the Newton method for the state system.
By default this is set to newton_verbose = False
so that there is not too much noise in the terminal.
The upcoming parameters are used to define the behavior of a Picard iteration, that may be used if we have multiple variables.
Note
For a detailed discussion of how to use the Picard iteration to solve a coupled state system, we refer to Coupled Problems  Picard Iteration. Note, that this demo is written for optimal control problems, but the definition of the state system can be transferred analogously to shape optimization problems, too.
First, we have a boolean flag, set via
picard_iteration = False
which determines whether the Picard iteration is enabled or not. This defaults
to picard_iteration = False
, so that the Picard solver is disabled by default.
The following two parameters determine, analogously to above, the tolerances for the
Picard iteration
picard_rtol = 1e10
picard_atol = 1e12
The default values for these parameters are picard_rtol = 1e10
and
picard_atol = 1e12
. Moreover, note that the tolerances of the Newton solver are adjusted automatically in case
a Picard iteration is performedm, so that an inexact Picard iteration is used.
The maximum amout of iterations for the Picard iteration are set with
picard_iter = 10
The default value for this is given by picard_iter = 50
.
Finally, we can enable verbose output of the Picard iteration with the following boolean flag
picard_verbose = False
which is set to picard_verbose = False
by default.
Section OptimizationRoutine#
The section OptimizationRoutine also closely resembles the one for optimal control problems. Again, we will take a brief look at all parameters here
The first parameter that can be controlled via the config file is algorithm
, which is
set via
algorithm = lbfgs
There are three possible choices for this parameter for shape optimization problems, namely
algorithm = gd
oralgorithm = gradient_descent
: A gradient descent methodalgorithm = cg
,algorithm = conjugate_gradient
,algorithm = ncg
,algorithm = nonlinear_cg
: Nonlinear CG methodsalgorithm = lbfgs
oralgorithm = bfgs
: limited memory BFGS method.
Thereafter, we specify the tolerances for the optimization algorithm with the parameters
rtol = 5e3
atol = 0.0
Again, rtol
denotes the relative, and atol
the absolute tolerance, and the
defaults for these parameters are given by rtol = 1e3
, and atol = 0.0
.
The next parameter is used to control the maximum number of iterations performed by the optimization algorithm. It is set via
max_iter = 50
and defaults to max_iter = 100
.
Next up, we have the initial guess for the step size, which can be determined via
initial_stepsize = 1.0
The default behavior is given by initial_stepsize = 1.0
.
The next parameter is given by
safeguard_stepsize = True
This parameter can be used to activate safeguarding of the initial stepsize for line search methods. This helps to choose an apropriate stepsize for the initial iteration even if the problem is poorly scaled.
The upcoming parameters are used for the Armijo rule
epsilon_armijo = 1e4
beta_armijo = 2
They are used to verify that the condition
holds, and if this is not satisfied, the stepsize is updated via \(t = \frac{t}{\beta_{\text{Armijo}}}\).
As default values for these parameters we use epsilon_armijo = 1e4
as well
as beta_armijo = 2
.
Next, we have a set of two parameters which detail the methods used for computing gradients in cashocs. These parameters are
gradient_method = direct
as well as
gradient_tol = 1e9
The first parameter, gradient_method
can be either gradient_method = direct
or gradient_method = iterative
. In the former case, a
direct solver is used to compute the gradient (using a Riesz projection) and in the latter case, an
iterative solver is used to do so. In case we have gradient_method = iterative
, the parameter
gradient_tol
is used to specify the (relative) tolerance for the iterative solver, in the other case
the parameter is not used.
The following parameter, soft_exit
, is a boolean flag which determines how
the optimization algorithm is terminated in case it does not converge. If soft_exit = True
, then an
error message is printed, but code after the solve
call of the
optimization problem will still be executed. However, when soft_exit = False
, cashocs
raises an exception and terminates. This is set via
soft_exit = False
and is set to soft_exit = False
by default.
Section LineSearch#
In this section, parameters regarding the line search can be specified. The type of the line search can be chosen via the parameter
method = armijo
Possible options are method = armijo
, which performs a simple backtracking line search based on the armijo rule with fixed steps (think of halving the stepsize in each iteration), and method = polynomial
, which uses polynomial models of the cost functional restricted to the line to generate “better” guesses for the stepsize. The default is method = armijo
.
The next parameter, polynomial_model
, specifies, which type of polynomials are used to generate new trial stepsizes. It is set via
polynomial_model = cubic
The parameter can either be polynomial_model = quadratic
or polynomial_model = cubic
. If this is polynomial_model = quadratic
, a quadratic interpolation polynomial along the search direction is generated and this is minimized analytically to generate a new trial stepsize. Here, only the current function value, the direction derivative of the cost functional in direction of the search direction, and the most recent trial stepsize are used to generate the polynomial. In case that polynomial_model = cubic
, the last two trial stepsizes (when available) are used in addition to the current cost functional value and the directional derivative, to generate a cubic model of the onedimensional cost functional, which is then minimized to compute a new trial stepsize.
For the polynomial models, we also have a safeguarding procedure, which ensures that trial stepsizes cannot be chosen too large or too small, and which can be configured with the following two parameters. The trial stepsizes generate by the polynomial models are projected to the interval \([\beta_{low} \alpha, \beta_{high} \alpha]\), where \(\alpha\) is the previous trial stepsize and \(\beta_{low}, \beta_{high}\) are factors which can be set via the parameters factor_low
and factor_high
. In the config file, this can look like this
factor_high = 0.5
factor_low = 0.1
and the values specified here are also the default values for these parameters.
Finally, we have the parameter
fail_if_not_converged = False
which determines, whether the line search is terminated if the state system cannot be solved at the current iterate. If this is fail_if_not_converged = True
, then an exception is raised. Otherwise, the iterate is counted as having too high of a function value and the stepsize is “halved” and a new iterate is formed.
Section AlgoLBFGS#
Next, we discuss the parameters relevant for the limited memory BFGS method. For details regarding this method, we refer to Schulz, Siebenborn, and Welker, Efficient PDE Constrained Shape Optimization Based on SteklovPoincaréType Metrics, where the methods are introduced.
The first parameter, bfgs_memory_size
, determines how large the storage of the BFGS method is. It is set via
bfgs_memory_size = 3
Usually, a higher storage leads to a better Hessian approximation, and thus to faster
convergence. However, this also leads to an increased memory usage. Typically, values
below 5 already work very well. The default is bfgs_memory_size = 5
.
The other parameter for the BFGS method is
use_bfgs_scaling = True
This determines, whether one should use a scaling of the initial Hessian approximation (see Nocedal and Wright, Numerical Optimization). This is usually very beneficial and should be kept enabled (which is the default).
Third, we have the parameter bfgs_periodic_restart
, which is set in the line
bfgs_periodic_restart = 0
This is a nonnegative integer value, which indicates the number of BFGS iterations, before a reinitialization takes place. In case that this is bfgs_periodic_restart = 0
(which is the default), no restarts are performed.
Finally, we have the parameter damped
, which can be set with
damped = False
This parameter is a boolean flag, which indicates whether Powell’s damping (on H) should be used or not. This is useful, when the curvature condition is not satisfied and (without damping) a restart would be required. The default is damped = False
.
Section AlgoCG#
The following parameters are used to define the behavior of the nonlinear conjugate gradient methods for shape optimization. For more details on this, we refer to the preprint Blauth, Nonlinear Conjugate Gradient Methods for PDE Constrained Shape Optimization Based on SteklovPoincaréType Metrics.
First, we define which nonlinear CG method is used by
cg_method = DY
Available options are
cg_method = FR
: The FletcherReeves methodcg_method = PR
: The PolakRibiere methodcg_method = HS
: The HestenesStiefel methodcg_method = DY
: The DaiYuan methodcg_method = HZ
: The HagerZhang method
The default value is cg_method = FR
. As for optimal control problems, the subsequent parameters are used to define the
restart behavior of the nonlinear CG methods. First, we have
cg_periodic_restart = False
This boolean flag en or disables that the NCG methods are restarted after a fixed amount of iterations, which is specified via
cg_periodic_its = 5
i.e., if cg_periodic_restart = True
and cg_periodic_its = n
, then the NCG method
is restarted every \(n\) iterations. The default behavior is given by
cg_periodic_restart = False
and cg_periodic_its = 10
.
Alternatively, there also exists a relative restart criterion (see Nocedal and Wright,
Numerical Optimization), which can be enabled
via the boolean flag cg_relative_restart
, which is defined in the line
cg_relative_restart = False
and the corresponding restart tolerance is set in
cg_restart_tol = 0.5
Note, that cg_restart_tol
should be in \((0, 1)\). If two subsequent
gradients generated by the nonlinear CG method are not “sufficiently
orthogonal”, the method is restarted with a gradient step. The default behavior
is given by cg_relative_restart = False
and cg_restart_tol = 0.25
.
Section ShapeGradient#
After we have specified the behavior of the solution algorithm, this section is used to specify parameters relevant to the computation of the shape gradient. Note, that by shape gradient we refer to the following object.
Let \(\mathcal{S} \subset \{ \Omega \;\vert\; \Omega \subset \mathbb{R}^d \}\) be a subset of the power set of \(\mathbb{R}^d\). Let \(J\) be a shape differentiable functional \(J \colon \mathcal{S} \to \mathbb{R}\) with shape derivative \(dJ(\Omega)[\mathcal{V}]\). Moreover, let \(a \colon H \times H \to \mathbb{R}\) be a symmetric, continuous, and coercive bilinear form on the Hilbert space \(H\). Then, the shape gradient \(\mathcal{G}\) of \(J\) (w.r.t. \(a\)) is defined as the solution of the problem
For PDE constrained shape optimization, it is common to use a bilinear form based on the linear elasticity equations, which enables smooth mesh deformations. This bilinear form is given as follows, in a general form, that is also implemented in cashocs
where \(H\) is some suitable subspace of \(H^1(\Omega)^d\) and \(\varepsilon(\mathcal{V}) = \frac{1}{2}(D\mathcal{V} + D\mathcal{V}^\top)\) is the symmetric part of the Jacobian. The subspace property is needed to include certain geometrical constraints of the shape optimization problem, which fix certain boundaries, into the shape gradient. For a detailed description of this setting we refer to the preprint Blauth, Nonlinear Conjugate Gradient Methods for PDE Constrained Shape Optimization Based on SteklovPoincaréType Metrics. Moreover, we note that for the second Lamé parameter \(\mu\), cashocs implements an idea from Schulz and Siebenborn, Computational Comparison of Surface Metric for PDE Constrained Shape Optimization: There, it is proposed to compute \(\mu\) as the solution of the Laplace problem
This allows to give the deformable and fixed boundaries a different stiffness, which is then smoothly extended into the interior of the domain. Moreover, they propose to use the solution of this Laplace equation directly for 2D problems, and to use \(\sqrt{\mu}\) for 3D problems.
Moreover, let us take a look at the possible _typing of boundaries that can be used with cashocs. In principle, there exist two _typing: deformable and fixed boundaries. On fixed boundaries, we impose homogeneous Dirichlet boundary conditions for the shape gradient, so that these are not moved under the corresponding deformation. In cashocs, we define what boundaries are fixed and deformable via their markers, which are either defined in the corresponding python script, or in the GMSH file, if such a mesh is imported.
The config file for Shape Optimization with a Poisson Problem defines the deformable boundaries with the command
shape_bdry_def = [1]
Note
Remember, that in Shape Optimization with a Poisson Problem, we defined boundaries
with the commands
boundary = CompiledSubDomain('on_boundary')
boundaries = MeshFunction('size_t', mesh, dim=1)
boundary.mark(boundaries, 1)
Hence, we see that the marker 1
corresponds to the entire boundary, so that this
is set to being deformable through the config.
As we do not have a fixed boundary for this problem, the corresponding list for the fixed boundaries is empty
shape_bdry_fix = []
Note, that cashocs also gives you the possibility of defining partially constrainted boundaries, where only one axial component is fixed, whereas the other two are not. These are defined in
shape_bdry_fix_x = []
shape_bdry_fix_y = []
shape_bdry_fix_z = []
For these, we have that shape_bdry_fix_x
is a list of all markers whose corresponding
boundaries should not be deformable in xdirection, but can be deformed in the y
and zdirections. Of course you can constrain a boundary to be only variable in a
single direction by adding the markers to the remaining lists.
The next parameter is specified via
use_pull_back = True
This parameter is used to determine, whether the material derivative should be computed for objects that are not state or adjoint variables. This is enabled by default.
Warning
This parameter should always be set to use_pull_back = True
, otherwise the shape derivative might
be wrong. Only disable it when you are sure what you are doing.
Furthermore, note that the material derivative computation is only correct, as long as no differential operators act on objects that are not state or adjoint variables. However, this should only be a minor restriction and not relevant for almost all problems.
Note
See Inverse Problem in Electric Impedance Tomography for a case, where we use
use_pull_back = False
.
The next parameters determine the coefficients of the bilinear form \(a\). First, we have the first Lamé parameter \(\lambda\), which is set via
lambda_lame = 1.428571428571429
The default value for this is lambda_lame = 0.0
.
Next, we specify the damping parameter \(\delta\) with the line
damping_factor = 0.2
The default for this is damping_factor = 0.0
.
Note
As the default value for the damping factor is damping_factor = 0.0
, this
should be set to a positive value in case the entire boundary of a problem
is deformable. Otherwise, the Riesz identification problem for the shape
gradient is not wellposed.
Finally, we define the values for \(\mu_\text{def}\) and \(\mu_\text{fix}\) via
mu_fix = 0.35714285714285715
mu_def = 0.35714285714285715
The default behavior is given by mu_fix = 1.0
and mu_def = 1.0
.
The parameter use_sqrt_mu
is a boolean flag, which switches between using
\(\mu\) and \(\sqrt{\mu}\) as the stiffness for the linear elasticity
equations, as discussed above. This is set via
use_sqrt_mu = False
and the default value is use_sqrt_mu = False
.
The next line in the config file is
inhomogeneous = False
This determines, whether an inhomogeneous linear elasticity equation is used to project the shape gradient. This scales the parameters \(\mu, \lambda\) and \(\delta\) by \(\frac{1}{\text{vol}}\), where \(\text{vol}\) is the volume of the current element (during assembly). This means, that smaller elements get a higher stiffness, so that the deformation takes place in the larger elements, which can handle larger deformations without reducing their quality too much. For more details on this approach, we refer to the paper Blauth, Leithäuser, and Pinnau, Model Hierarchy for the Shape Optimization of a Microchannel Cooling System.
Moreover, the parameter
update_inhomogeneous = False
can be used to update the local mesh size after each mesh deformation, in case this is update_inhomogeneous = True
, so that elements which become smaller also obtain a higher stiffness and vice versa. The default is update_inhomogeneous = False
.
For the inhomogeneous mesh stiffness, we also have the parameter inhomogeneous_exponent
, which is specified via
inhomogeneous_exponent = 1.0
This parameter can be used to specify an exponent for the inhomogeneous mesh stiffness, so that the parameters
\(\mu, \lambda\) and \(\delta\) are scaled by \(\left( \frac{1}{\text{vol}} \right)^p\), where
\(p\) is specified in inhomogeneous_exponent
. The default for this parameter is inhomogeneous_exponent = 1.0
.
There is also a different possibility to define the stiffness parameter \(\mu\) using cashocs, namely to define \(\mu\) in terms of how close a point of the computational domain is to a boundary. In the following we will explain this alternative way of defining \(\mu\). To do so, we must first set the boolean parameter
use_distance_mu = True
which enables this formulation and deactivates the previous one. Note that by default,
the value is use_distance_mu = False
. Next, we have the parameters dist_min
, dist_max
,
mu_min
and mu_max
. These do the following: If the distance to the boundary is
smaller than dist_min
, the value of \(\mu\) is set to mu_min
, and if the distance
to the boundary is larger than dist_max
, \(\mu\) is set to mu_max
. If the distance
to the boundary is between dist_min
and dist_max
, the value of \(\mu\) is
interpolated between mu_min
and mu_max
. The type of this interpolation is
determined by the parameter
smooth_mu = True
If this parameter is set to smooth_mu = True
, then a smooth, cubic polynomial is used to
interplate between mu_min
and mu_max
, which yields a continuously differentiable
\(\mu\). If this is set to smooth_mu = False
, then a linear interpolation is used, which only yields
a continuous \(\mu\). The default for this parameter is smooth_mu = False
.
Finally, we can specify which boundaries we want to incorporate when computing the distance. To do so, we can specify a list of indices which contain the boundary markers in the parameter
boundaries_dist = [1,2,3]
This means, that only boundaries marked with 1, 2, and 3 are considered for computing the distance, and all others are ignored. The default behavior is that all (outer) boundaries are considered.
There is also another possibility to compute the shape gradient in cashocs, namely using the \(p\)Laplacian, as proposed by Müller, Kühl, Siebenborn, Deckelnick, Hinze, and Rung. In order to do so, we have the following line
use_p_laplacian = False
If this is set to use_p_laplacian = True
, the \(p\)Laplacian is used to compute the shape gradient, as explained in Shape Optimization with the pLaplacian. However, by default this is disabled.
The value of \(p\) which is then used is defined in the next line
p_laplacian_power = 6
which defaults to p_laplacian_power = 2
, whenever the parameter is not defined. The higher \(p\) is chosen, the better the numerical are expected to be, but the numerical solution of the problem becomes more involved.
Finally, there is the possibility to use a stabilized weak form for the \(p\)Laplacian operator, where the stabilization parameter can be defined in the line
p_laplacian_stabilization = 0.0
The default value of this parameter is p_laplacian_stabilization = 0.0
. Note, that the parameter should be chosen comparatively small, i.e., significantly smaller than 1.
Furthermore, we have the parameter fixed_dimensions
, which enables us to restrict the shape gradient to specific dimensions. It is set via
fixed_dimensions = []
In case fixed_dimensions = []
, there is no restriction on the shape gradient. However, if fixed_dimensions = [i]
, then the ith component of the shape gradient is set to 0, so that we have no deformation in the ith coordinate direction. For example, if fixed_dimensions = [0, 2]
, we only have a deformation in the \(y\)component of the mesh. The default is fixed_dimensions = []
.
Section Regularization#
In this section, the parameters for shape regularizations are specified. For a detailed discussion of their usage, we refer to Regularization for Shape Optimization Problems.
First, we have the parameters factor_volume
and target_volume
. These are set
via the lines
factor_volume = 0.0
target_volume = 3.14
They are used to implement the (target) volume regularization term
Here, \(\mu_\text{vol}\) is specified via factor_volume
, and \(\text{vol}_\text{des}\)
is the target volume, specified via target_volume
. The default behavior is
factor_volume = 0.0
and target_volume = 0.0
, so that we do not have
a volume regularization.
The next line, i.e.,
use_initial_volume = True
determines the boolean flag use_initial_volume
. If this is set to use_initial_volume = True
,
then not the value given in target_volume
is used, but instead the
volume of the initial geometry is used for \(\text{vol}_\text{des}\).
For the next two _typing of regularization, namely the (target) surface and (target) barycenter regularization, the syntax for specifying the parameters is completely analogous. For the (target) surface regularization we have
factor_surface = 0.0
target_surface = 1.0
These parameter are used to implement the regularization term
Here, \(\mu_\text{surf}\) is determined via factor_surface
, and
\(\text{surf}_\text{des}\) is determined via target_surface
. The default
values are given by factor_surface = 0.0
and target_surface = 0.0
.
As for the volume regularization, the parameter
use_initial_surface = True
determines whether the target surface area is specified via target_surface
or if the surface area of the initial geometry should be used instead. The default
behavior is given by use_initial_surface = False
.
Next, we have the curvature regularization, which is controlled by the parameter
factor_curvature = 0.0
This is used to determine the size of \(\mu_\text{curv}\) in the regularization term
where \(\kappa\) denotes the mean curvature. This regularization term can be used to generate more smooth boundaries and to prevent kinks from occurring.
Finally, we have the (target) barycenter regularization. This is specified via the parameters
factor_barycenter = 0.0
target_barycenter = [0.0, 0.0, 0.0]
and implements the term
The default behavior is given by factor_barycenter = 0.0
and target_barycenter = [0,0,0]
,
so that we do not have a barycenter regularization.
The flag
use_initial_barycenter = True
again determines, whether \(\text{bary}_\text{des}\) is determined via target_barycenter
or if the barycenter of the initial geometry should be used instead. The default behavior
is given by use_initial_barycenter = False
.
Hint
The object target_barycenter
has to be a list. For 2D problems it is also
sufficient, if the list only has two entries, for the \(x\) and \(y\)
barycenters.
Finally, we have the parameter use_relative_scaling
which is set in the line
use_relative_scaling = False
This boolean flag does the following. For some regularization term \(J_\text{reg}(\Omega)\) with corresponding
factor \(\mu\) (as defined above), the default behavior is given by use_relative_scaling = False
adds the term \(\mu J_\text{reg}(\Omega)\) to the cost functional, so that the
factor specified in the configuration file is actually used as the factor for the regularization term.
In case use_relative_scaling = True
, the behavior is different, and the following term is
added to the cost functional: \(\frac{\mu}{\left\lvert J_\text{reg}(\Omega_0) \right\rvert} J_\text{reg}(\Omega)\),
where \(\Omega_0\) is the initial guess for the geometry. In particular, this means
that the magnitude of the regularization term is equal to \(\mu\) on the initial geometry.
This allows a detailed weighting of multiple regularization terms, which is particularly
useful in case the cost functional is also scaled (see Scaling of the Cost Functional).
Section MeshQuality#
This section details the parameters that influence the quality of the computational mesh. First, we have the lines
volume_change = inf
angle_change = inf
These parameters are used to specify how much the volume and the angles, respectively, of the mesh elements are allowed to change in a single transformation. In particular, they implement the following criteria (see Etling, Herzog, Loayza, Wachsmuth, First and Second Order Shape Optimization Based on Restricted Mesh Deformations)
Here, \(\alpha\) corresponds to volume_change
and \(\beta\) corresponds
to angle_change
, and \(\mathcal{V}\) is the deformation. The default behavior
is given by volume_change = inf
and angle_change = inf
, so that no restrictions
are posed. Note, that, e.g., Etling, Herzog, Loayza, Wachsmuth,
First and Second Order Shape Optimization Based on Restricted Mesh Deformations use the values volume_change = 2.0
and
angle_change = 0.3
.
The next two parameters are given byx
tol_lower = 0.0
tol_upper = 1e15
These parameters specify a kind of interval for the mesh quality. In particular, we have the following situation (note that the mesh quality is always an element in \([0,1]\)):
If the mesh quality is in \([\texttt{tol upper}, 1]\), the mesh is assumed to be “good”, so that finite element solutions of the corresponding PDEs are sensible and not influenced by the mesh quality or discretization artifacts.
If the mesh quality is in \([\texttt{tol lower}, \texttt{tol upper}]\), a kind of breaking point is reached. Here, it is assumed that the mesh is sufficiently good so that the solution of the state system is still possible. However, a mesh whose quality is in this interval should not be used anymore to compute the solution of the adjoint system or to compute the shape gradient, as the quality is too poor for this purpose. Usually, this means that the algorithm is terminated, unless remeshing is enabled. In the latter case, remeshing is performed.
If the mesh quality is in the interval \([0, \texttt{tol lower}]\), the mesh quality is assumed to be so poor, that even the solution of the state system is not possible anymore. In practice, this can only happen during the Armijo line search. Thanks to our previous considerations, we also know that the mesh, that is to be deformed, has at least a quality of
tol_lupper
, so that this quality might be reached again, if the step size is just decreased sufficiently often. This way, it is ensured that the state system is only solved when the mesh quality is larger thantol_lower
, so that the corresponding cost functional value is reasonable.
The default behavior is given by tol_lower = 0.0
and tol_upper = 1e15
,
so that there are basically no requirements on the mesh quality.
Finally, the upcoming two parameters specify how exactly the mesh quality is measured. The first one is
measure = condition_number
and determines one of the four mesh quality criteria, as defined in MeshQuality
.
Available options are
measure = skewness
measure = maximum_angle
measure = radius_ratios
measure = condition_number
(see MeshQuality
for a detailed description).
The default value is given by measure = skewness
.
The parameter type
determines, whether the minimum quality over all
elements (type = min
) or the average quality over all elements (type = avg
)
shall be used. This is set via
type = min
and defaults to type = min
.
Finally, we have the parameter remesh_iter
in which the user can specify after how many iterations a remeshing should be performed. It is given by
remesh_iter = 0
where remesh_iter = 0
means that no automatic remeshing is performed (this is the default), and remesh_iter = n
means that remeshing is performed after each n iterations. Note that to use this parameter and avoid unexpected results, it might be beneficial to the the lower and upper mesh quality tolerances to a low value, so that the “quality based remeshing” does not interfere with the “iteration based remeshing”, but both can be used in combination.
Section Output#
In this section, the parameters for the output of the algorithm, either in the terminal
or as files, are specified. First, we have the parameter verbose
. This is used to toggle the output of the
optimization algorithm. It defaults to verbose = True
and is controlled via
verbose = True
The parameter save_results
is a boolean flag, which determines whether a history
of the optimization algorithm, including cost functional value, gradient norm, accepted
step sizes, and mesh quality, shall be saved to a .json file. This defaults to save_results = True
,
and can be set with
save_results = False
Moreover, we define the parameter save_txt
save_txt = False
This saves the output of the optimization, which is usually shown in the terminal, to a .txt file, which is humanreadable.
The next line in the config file is
save_state = False
Here, the parameter save_state
is set. This is a boolean flag, which can be set to
save_state = True
to enable that cashocs generates .xdmf files for the state variables for each iteration the optimization algorithm performs. These are great for visualizing the
steps done by the optimization algorithm, but also need some disc space, so that they are disabled by default.
Note, that for visualizing these files, you need Paraview.
The next parameter, save_adjoint
works analogously, and is given in the line
save_adjoint = False
If this is set to True, cashocs generates .xdmf files for the adjoint variables in each iteration of the optimization algorithm. Its main purpose is for debugging.
The next parameter is given by save_gradient
, which is given in the line
save_gradient = False
This boolean flag ensures that a paraview with the computed shape gradient is saved in result_dir/xdmf
. The main purpose of this is for debugging.
Moreover, we also have the parameter save_mesh
that is set via
save_mesh = False
This is used to save the optimized geometry to a GMSH file. The default behavior
is given by save_mesh = False
. Note, that this is only
possible if the input mesh was already generated by GMSH, and specified in the Mesh
section of the config file. For any other meshes, the underlying mesh is also saved in
the .xdmf files, so that you can at least always visualize the optimized geometry.
In the end, we also have, like for optimal control problems, a parameter that specifies
where the output is placed, again named result_dir
, which is given in the config file
in the line
result_dir = ./results
As before, this is either a relative or absolute path to the directory where the results should be placed.
The parameter precision
, which is set via
precision = 3
is an integer parameter which determines how many significant digits are printed in the output to the console and / or the result file.
Moreover, we have the parameter time_suffix
, which adds a suffix to the result directory based on the current time. It is controlled by the line
time_suffix = False
Summary#
Finally, an overview over all parameters and their default values can be found in the following.
[Mesh]#
Parameter = Default value 
Remarks 


Only needed for remeshing 

Only needed for remeshing 

Only needed for remeshing 

if 

if 
[StateSystem]#
Parameter = Default value 
Remarks 


using 

relative tolerance for Newton’s method 

absolute tolerance for Newton’s method 

maximum iterations for Newton’s method 

if 

if 





relative tolerance for Picard iteration 

absolute tolerance for Picard iteration 

maximum iterations for Picard iteration 


[OptimizationRoutine]#
Parameter = Default value 
Remarks 


has to be specified by the user; see 

relative tolerance for the optimization algorithm 

absolute tolerance for the optimization algorithm 

maximum iterations for the optimization algorithm 

specifies the solver for computing the gradient, can be either 

the relative tolerance in case an iterative solver is used to compute the gradient. 

if 
[LineSearch]#
Parameter = Default value 
Remarks 




initial stepsize for the first iteration in the Armijo rule 





De(activates) a safeguard against poor scaling 

This specifies, whether a cubic or quadratic model is used for computing trial stepsizes 

Safeguard for stepsize, upper bound 

Safeguard for stepsize, lower bound 

if this is 
[AlgoLBFGS]#
Parameter = Default value 
Remarks 


memory size of the LBFGS method 

if 

specifies, after how many iterations the method is restarted. If this is 0, no restarting is done. 

specifies whether damping for the BFGS method should be used to enforce the curvature condition and prevent restarting 
[AlgoCG]#
Parameter = Default value 
Remarks 


specifies which nonlinear CG method is used 

if 

specifies, after how many iterations the NCG method is restarted, if applicable 

if 

the tolerance of the relative restart criterion, if applicable 
[ShapeGradient]#
Parameter = Default value 
Remarks 


list of indices for the deformable boundaries 

list of indices for the fixed boundaries 

list of indices for boundaries with fixed x values 

list of indices for boundaries with fixed y values 

list of indices for boundaries with fixed z values 

a list of coordinates which should be fixed during the shape optimization (x=0, y=1, etc.) 

if 

value of the first Lamé parameter for the elasticity equations 

value of the damping parameter for the elasticity equations 

value of the second Lamé parameter on the deformable boundaries 

value of the second Lamé parameter on the fixed boundaries 

if 

if 

if 

The exponent for the inhomogeneous mesh stiffness 

if 

Specifies the distance to the boundary, until which \(\mu\) is given by 

Specifies the distance to the boundary, until which \(\mu\) is given by 

The value of \(\mu\) for a boundary distance smaller than 

The value of \(\mu\) for a boundary distance larger than 

The indices of the boundaries, which shall be used to compute the distance, 

If false, a linear (continuous) interpolation between 

If 

The parameter \(p\) of the \(p\)Laplacian 

The stabilization parameter for the \(p\)Laplacian problem. No stabilization is used when this is 
[Regularization]#
Parameter = Default value 
Remarks 


value of the regularization parameter for volume regularization; needs to be nonnegative 

prescribed volume for the volume regularization 

if 

value of the regularization parameter for surface regularization; needs to be nonnegative 

prescribed surface for the surface regularization 

if 

value of the regularization parameter for curvature regularization; needs to be nonnegative 

value of the regularization parameter for barycenter regularization; needs to be nonnegative 

prescribed barycenter for the barycenter regularization 

if 

if 
[MeshQuality]#
Parameter = Default value 
Remarks 


determines by what factor the volume of a cell is allowed to change within a single deformation 

determines how much the angles of a cell are allowed to change within a single deformation 

if the mesh quality is lower than this tolerance, the state system is not solved for the Armijo rule, instead step size is decreased 

if the mesh quality is between 

determines which quality measure is used 

determines if minimal or average quality is considered 
[Output]#
Parameter = Default value 
Remarks 


if 

if 

if 

if 

if 

if 

if 

path to the directory, where the output should be placed 

number of significant digits to be printed 

Boolean flag, which adds a suffix to 