Stokes integrators
stokes_integrators
Solver classes targetting systems dealing with momentum conservation.
This module provides a fine-tuned abstract base class, StokesSolverBase
, from which
efficient numerical solvers can be derived, such as StokesSolver
for the Stokes system
of conservation equations and ViscoelasticStokesSolver
for an incremental displacement
formulation of the latter using a Maxwell rheology. The module also exposes a function,
create_stokes_nullspace
, to automatically generate null spaces compatible with PETSc
solvers and a class, BoundaryNormalStressSolver
, to solve for the normal stress acting
on a domain boundary. Typically, users instantiate the StokesSolver
class by providing
a relevant set of arguments and then call the solve
method to request a solver update.
iterative_stokes_solver_parameters = {'mat_type': 'matfree', 'ksp_type': 'preonly', 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'pc_fieldsplit_schur_type': 'full', 'fieldsplit_0': {'ksp_type': 'cg', 'ksp_rtol': 1e-05, 'ksp_max_it': 1000, 'pc_type': 'python', 'pc_python_type': 'gadopt.SPDAssembledPC', 'assembled_pc_type': 'gamg', 'assembled_mg_levels_pc_type': 'sor', 'assembled_pc_gamg_threshold': 0.01, 'assembled_pc_gamg_square_graph': 100, 'assembled_pc_gamg_coarse_eq_limit': 1000, 'assembled_pc_gamg_mis_k_minimum_degree_ordering': True}, 'fieldsplit_1': {'ksp_type': 'fgmres', 'ksp_rtol': 0.0001, 'ksp_max_it': 200, 'pc_type': 'python', 'pc_python_type': 'firedrake.MassInvPC', 'Mp_pc_type': 'ksp', 'Mp_ksp_ksp_rtol': 1e-05, 'Mp_ksp_ksp_type': 'cg', 'Mp_ksp_pc_type': 'sor'}}
module-attribute
Default iterative solver parameters for solution of Stokes system.
We configure the Schur complement approach as described in Section of 4.3 of Davies et al. (2022), using PETSc's fieldsplit preconditioner type, which provides a class of preconditioners for mixed problems that allows a user to apply different preconditioners to different blocks of the system.
The fieldsplit_0
entries configure solver options for the velocity block. The linear
systems associated with this matrix are solved using a combination of the Conjugate
Gradient (cg
) method and an algebraic multigrid preconditioner (gamg
).
The fieldsplit_1
entries contain solver options for the Schur complement solve itself.
For preconditioning, we approximate the Schur complement matrix with a mass matrix
scaled by viscosity, with the viscosity sourced from the approximation object provided
to Stokes solver. Since this preconditioner step involves an iterative solve, the Krylov
method used for the Schur complement needs to be of flexible type, and we use the
Flexible Generalized Minimal Residual (fgmres
) method by default.
Note
G-ADOPT will use the above default iterative solver parameters if the argument
solver_parameters="iterative"
is provided or in 3-D if the solver_parameters
argument is omitted. To make modifications to these default values, the most
convenient approach is to provide the modified values as a dictionary via the
solver_parameters_update
argument. This dictionary can also hold new pairs of keys
and values to extend the default ones.
.
direct_stokes_solver_parameters = {'mat_type': 'aij', 'ksp_type': 'preonly', 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps'}
module-attribute
Default direct solver parameters for solution of Stokes system.
We configure the direct solver to use the LU (lu
) factorisation provided by the MUMPS
package.
Note
G-ADOPT will use the above default direct solver parameters if the argument
solver_parameters="direct"
is provided or in 2-D if the solver_parameters
argument is omitted. To make modifications to these default values, the most
convenient approach is to provide the modified values as a dictionary via the
solver_parameters_update
argument. This dictionary can also hold new pairs of keys
and values to extend the default ones.
newton_stokes_solver_parameters = {'snes_type': 'newtonls', 'snes_linesearch_type': 'l2', 'snes_max_it': 100, 'snes_atol': 1e-10, 'snes_rtol': 1e-05}
module-attribute
Default non-linear solver parameters for solution of Stokes system.
We use a setup based on Newton's method (newtonls) with a secant line search over the L2-norm of the function.
Note
G-ADOPT will use the above default non-linear solver parameters in conjunction with
the above iterative or default solver parameters if the viscosity, mu
, provided to
the approximation depends on the solver's solution. To make modifications to these
default values, the most convenient approach is to provide the modified values as a
dictionary via the solver_parameters_update
argument. This dictionary can also hold
new pairs of keys and values to extend the default ones.
MetaPostInit
Bases: ABCMeta
Calls the implemented prepare_solver
method after __init__
returns.
The implemented behaviour allows any subclass __init__
method to first call its
parent class's __init__
through super(), then execute its own code, and finally
call prepare_solver
. The latter call is automatic and does not require any
attention from the developer or user.
StokesSolverBase(solution, approximation, /, *, dt=None, theta=0.5, additional_forcing_term=None, bcs={}, quad_degree=6, solver_parameters=None, solver_parameters_update=None, J=None, constant_jacobian=False, nullspace=None, transpose_nullspace=None, near_nullspace=None)
Bases: ABC
Solver for a system involving mass and momentum conservation.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
solution
|
Function
|
Firedrake function representing the field over the mixed Stokes space |
required |
approximation
|
BaseApproximation
|
G-ADOPT approximation defining terms in the system of equations |
required |
dt
|
float | None
|
Float specifying the time step if the system involves a coupled time integration |
None
|
theta
|
float
|
Float defining the theta scheme parameter used in a coupled time integration |
0.5
|
additional_forcing_term
|
Form | None
|
Firedrake form specifying an additional term contributing to the residual |
None
|
bcs
|
dict[int | str, dict[str, Any]]
|
Dictionary specifying boundary conditions (identifier, type, and value) |
{}
|
quad_degree
|
int
|
Integer denoting the quadrature degree |
6
|
solver_parameters
|
dict[str, str | float] | str | None
|
Dictionary of PETSc solver options or string matching one of the default sets |
None
|
solver_parameters_update
|
dict[str, str | float] | str | None
|
Dictionary of PETSc solver options used to update the default G-ADOPT options |
None
|
J
|
Function | None
|
Firedrake function representing the Jacobian of the mixed Stokes system |
None
|
constant_jacobian
|
bool
|
Boolean specifying whether the Jacobian of the system is constant |
False
|
nullspace
|
MixedVectorSpaceBasis
|
A |
None
|
transpose_nullspace
|
MixedVectorSpaceBasis
|
A |
None
|
near_nullspace
|
MixedVectorSpaceBasis
|
A |
None
|
Valid keys for boundary conditions
Condition | Type | Description |
---|---|---|
u | Strong | Solution |
ux | Strong | Solution along the first Cartesian axis |
uy | Strong | Solution along the second Cartesian axis |
uz | Strong | Solution along the third Cartesian axis |
un | Weak | Solution along the normal to the boundary |
stress | Weak | Traction across the boundary |
normal_stress | Weak | Stress component normal to the boundary |
free_surface | Weak | Free-surface characteristics of the boundary |
Valid keys describing the free surface boundary
Argument | Required | Description |
---|---|---|
delta_rho_fs | Yes (d) | Density contrast along the free surface |
RaFS | Yes (nd) | Rayleigh number (free-surface density contrast) |
variable_rho_fs | No | Account for buoyancy effects on interior density |
Classic theta values for coupled implicit time integration
Theta | Scheme |
---|---|
0.5 | Crank-Nicolson method |
1.0 | Backward Euler method |
Note: Such a coupling can arise in the case of a free-surface implementation.
Source code in g-adopt/gadopt/stokes_integrators.py
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 |
|
prepare_solver()
Runs methods to set up the variational problem and solver.
Source code in g-adopt/gadopt/stokes_integrators.py
364 365 366 367 368 369 370 |
|
set_boundary_conditions()
Sets strong and weak boundary conditions.
Source code in g-adopt/gadopt/stokes_integrators.py
372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 |
|
set_free_surface_boundary(params_fs, bc_id)
Sets the given boundary as a free surface.
This method calculates the normal stress at the free surface boundary. In the
coupled approach, it also sets a zero-interior strong condition away from that
boundary and populates the free_surface_map
dictionary used to calculate the
free-surface contribution to the momentum weak form.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
params_fs
|
dict[str, int | bool]
|
Dictionary holding information about the free surface boundary |
required |
bc_id
|
int | str
|
Integer representing the index of the mesh boundary |
required |
Returns:
Type | Description |
---|---|
Expr
|
UFL expression for the normal stress at the free surface boundary |
Source code in g-adopt/gadopt/stokes_integrators.py
406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 |
|
set_equations()
abstractmethod
Sets Equation instances for each equation in the system.
Equations must be ordered like solutions in the mixed space.
Source code in g-adopt/gadopt/stokes_integrators.py
426 427 428 429 430 431 |
|
set_form()
Sets the weak form including linear and bilinear terms.
Source code in g-adopt/gadopt/stokes_integrators.py
433 434 435 436 437 438 439 440 441 442 443 444 |
|
set_petsc_options()
Sets PETSc solver options.
Source code in g-adopt/gadopt/stokes_integrators.py
446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 |
|
set_solver()
Sets up the Firedrake variational problem and solver.
Source code in g-adopt/gadopt/stokes_integrators.py
490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 |
|
solve()
Solves the system.
Source code in g-adopt/gadopt/stokes_integrators.py
526 527 528 529 |
|
StokesSolver(solution, approximation, T=0.0, /, **kwargs)
Bases: StokesSolverBase
Solver for the Stokes system.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
solution
|
Function
|
Firedrake function representing the field over the mixed Stokes space |
required |
approximation
|
BaseApproximation
|
G-ADOPT approximation defining terms in the system of equations |
required |
T
|
Function | float
|
Firedrake function representing the temperature field |
0.0
|
dt
|
Float quantifying the time step used in a coupled time integration |
required | |
theta
|
Float quantifying the implicit contribution in a coupled time integration |
required | |
additional_forcing_term
|
Firedrake form specifying an additional term contributing to the residual |
required | |
bcs
|
Dictionary specifying boundary conditions (identifier, type, and value) |
required | |
quad_degree
|
Integer denoting the quadrature degree |
required | |
solver_parameters
|
Dictionary of PETSc solver options or string matching one of the default sets |
required | |
solver_parameters_update
|
Dictionary of PETSc solver options used to update the default G-ADOPT options |
required | |
J
|
Firedrake function representing the Jacobian of the mixed Stokes system |
required | |
constant_jacobian
|
Boolean specifying whether the Jacobian of the system is constant |
required | |
nullspace
|
A |
required | |
transpose_nullspace
|
A |
required | |
near_nullspace
|
A |
required |
Source code in g-adopt/gadopt/stokes_integrators.py
571 572 573 574 575 576 577 578 579 580 581 582 583 584 |
|
force_on_boundary(subdomain_id, **kwargs)
Computes the force acting on a boundary.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
subdomain_id
|
int | str
|
The subdomain ID of a physical boundary. |
required |
Returns:
Type | Description |
---|---|
Function
|
The force acting on the boundary. |
Source code in g-adopt/gadopt/stokes_integrators.py
665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 |
|
ViscoelasticStokesSolver(solution, approximation, stress_old, displacement, **kwargs)
Bases: StokesSolverBase
Solves the Stokes system assuming a Maxwell viscoelastic rheology.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
solution
|
Function
|
Firedrake function representing the field over the mixed Stokes space |
required |
approximation
|
SmallDisplacementViscoelasticApproximation
|
G-ADOPT approximation defining terms in the system of equations |
required |
stress_old
|
Function
|
Firedrake function representing the deviatoric stress at the previous time step |
required |
displacement
|
Function
|
Firedrake function representing the total displacement |
required |
dt
|
Float quantifying the time step used in a coupled time integration |
required | |
theta
|
Float quantifying the implicit contribution in a coupled time integration |
required | |
additional_forcing_term
|
Firedrake form specifying an additional term contributing to the residual |
required | |
bcs
|
Dictionary specifying boundary conditions (identifier, type, and value) |
required | |
quad_degree
|
Integer denoting the quadrature degree |
required | |
solver_parameters
|
Dictionary of PETSc solver options or string matching one of the default sets |
required | |
solver_parameters_update
|
Dictionary of PETSc solver options used to update the default G-ADOPT options |
required | |
J
|
Firedrake function representing the Jacobian of the mixed Stokes system |
required | |
constant_jacobian
|
Boolean specifying whether the Jacobian of the system is constant |
required | |
nullspace
|
A |
required | |
transpose_nullspace
|
A |
required | |
near_nullspace
|
A |
required |
Source code in g-adopt/gadopt/stokes_integrators.py
808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 |
|
set_equations()
Sets up UFL forms for the viscoelastic Stokes equations residual.
Source code in g-adopt/gadopt/stokes_integrators.py
840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 |
|
BoundaryNormalStressSolver(stokes_solver, subdomain_id, **kwargs)
A class for calculating surface forces acting on a boundary.
This solver computes topography on boundaries using the equation:
where \(sigma_(rr)\) is defined as:
Instead of assuming a unit normal vector \(hat n\), this solver uses FacetNormal
from Firedrake to accurately determine the normal vectors, which is particularly
useful for complex meshes like icosahedron meshes in spherical simulations.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
stokes_solver
|
StokesSolver
|
The Stokes solver providing the necessary fields for calculating stress. |
required |
subdomain_id
|
str | int
|
The subdomain ID of a physical boundary. |
required |
**kwargs
|
Optional keyword arguments. You can provide: - solver_parameters (dict[str, str | float]): Parameters to control the variational solver. If not provided, defaults are chosen based on whether the Stokes solver is direct or iterative. |
{}
|
Source code in g-adopt/gadopt/stokes_integrators.py
925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 |
|
solve()
Solves a linear system for the force and applies necessary boundary conditions.
Returns:
Type | Description |
---|---|
The modified force after solving the linear system and applying boundary conditions. |
Source code in g-adopt/gadopt/stokes_integrators.py
950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 |
|
create_stokes_nullspace(Z, closed=True, rotational=False, translations=None, ala_approximation=None, top_subdomain_id=None)
Create a null space for the mixed Stokes system.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
Z
|
WithGeometry
|
Firedrake mixed function space associated with the Stokes system |
required |
closed
|
bool
|
Whether to include a constant pressure null space |
True
|
rotational
|
bool
|
Whether to include all rotational modes |
False
|
translations
|
Optional[list[int]]
|
List of translations to include |
None
|
ala_approximation
|
Optional[AnelasticLiquidApproximation]
|
AnelasticLiquidApproximation for calculating (non-constant) right nullspace |
None
|
top_subdomain_id
|
Optional[str | int]
|
Boundary id of top surface. Required when providing ala_approximation. |
None
|
Returns:
Type | Description |
---|---|
MixedVectorSpaceBasis
|
A Firedrake mixed vector space basis incorporating the null space components |
Source code in g-adopt/gadopt/stokes_integrators.py
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 |
|
ala_right_nullspace(W, approximation, top_subdomain_id)
Compute pressure nullspace for Anelastic Liquid Approximation.
Parameters:
Name | Type | Description | Default |
---|---|---|---|
W
|
WithGeometry
|
pressure function space |
required |
approximation
|
AnelasticLiquidApproximation
|
AnelasticLiquidApproximation with equation parameters |
required |
top_subdomain_id
|
str | int
|
boundary id of top surface |
required |
Returns:
Type | Description |
---|---|
pressure nullspace solution |
To obtain the pressure nullspace solution for the Stokes equation in Anelastic Liquid Approximation, which includes a pressure-dependent buoyancy term, we try to solve the equation:
Taking the divergence:
then testing it with q:
followed by integration by parts:
This elliptic equation can be solved with natural boundary conditions by imposing our original equation above, which eliminates all boundary terms:
However, if we do so on all boundaries we end up with a system that has the same nullspace, as the one we are after (note that we ended up merely testing the original equation with \(nabla q\)). Instead we use the fact that the gradient of the null mode is always vertical, and thus the null mode is constant at any horizontal level (geoid), specifically the top surface. Choosing any nonzero constant for this surface fixes the arbitrary scalar multiplier of the null mode. We choose the value of one and apply it as a Dirichlet boundary condition.
Note that this procedure does not necessarily compute the exact nullspace of the discretised Stokes system. In particular, since not every test function \(v in V\), the velocity test space, can be written as \(v=nabla q\) with \(q in W\), the pressure test space, the two terms do not necessarily exactly cancel when tested with \(v\) instead of \(nabla q\) as in our final equation. However, in practice the discrete error appears to be small enough, and providing this nullspace gives an improved convergence of the iterative Stokes solver.
Source code in g-adopt/gadopt/stokes_integrators.py
686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 |
|