pyoomph.generic package

Submodules

Module contents

class pyoomph.generic.Equations(*args: Any, **kwargs: Any)[source]

Bases: BaseEquations

Equations to be solved on a domain, i.e. including spatial coordinates. Add unknown fields by overriding the define_fields() method and residuals by overriding the define_residuals() method.

See BaseEquations for further methods.

activate_coordinates_as_dofs(coordinate_space=None)[source]

Activates the coordinates as degrees of freedom (dofs) for a moving mesh. You must then add residuals in the define_residuals method for the field “mesh”, e.g.

def define_residuals(self):
x,xtest=var_and_test(“mesh”) X=var(“lagrangian”) self.add_weak(grad(x-X,lagrangian=True),grad(xtest,lagrangian=True),lagrangian=True)

for a Laplace-smoothed mesh

Parameters:: coordinate_space (Optional[str]) – The coordinate space to be set as the coordinate space for the element. Valid options are “C2TB”, “C2”, “C1TB”, or “C1”. If not provided, the coordinate space will not be set.
Raises:: ValueError – If the provided coordinate space is not one of the valid options.
Return type:: None
Returns:: None

define_scalar_field(name, space, scale=None, testscale=None, discontinuous_refinement_exponent=None)[source]

Define a scalar field on this domain. Must be called within the specified implementation of the method define_fields().

Parameters:

name (str) – The name of the vector field.
space (FiniteElementSpaceEnum) – The finite element space on which the vector field is defined.
scale (ExpressionNumOrNone) – The scale for the vector field for nondimensionalization. Defaults to None.
testscale (ExpressionNumOrNone) – The scale for the test function of the vector field for nondimensionalization. Defaults to None.
discontinuous_refinement_exponent (Optional[float]) – The exponent for the discontinuous refinement. Defaults to None.

define_vector_field(name, space, dim=None, scale=None, testscale=None)[source]

Define a vector field on this domain. Must be called within the specified implementation of the method define_fields().

Parameters:

name (str) – The name of the vector field.
space (FiniteElementSpaceEnum) – The finite element space on which the vector field is defined.
dim (Optional[int]) – The dimension of the vector field. If not provided, it defaults to the nodal dimension.
scale (ExpressionNumOrNone) – The scale for the vector field for nondimensionalization. Defaults to None.
testscale (ExpressionNumOrNone) – The scale for the test function of the vector field for nondimensionalization. Defaults to None.

get_opposite_parent_domain(raise_error_if_none=True)[source]

Returns the bulk domain at the opposite side of this interface.

Parameters:: raise_error_if_none (bool) – If there is no opposite side set, raise an error. Otherwise, just None is returned.
Return type:: Optional[FiniteElementCodeGenerator]
Returns:: The bulk domain at the opposite side.

get_opposite_side_of_interface(raise_error_if_none=True)[source]

Returns the interface domain at the opposite side of this interface.

Parameters:: raise_error_if_none (bool) – If there is no opposite side set, raise an error. Otherwise, just None is returned.
Return type:: Optional[FiniteElementCodeGenerator]
Returns:: The interface domain at the opposite side.

get_weak_dirichlet_terms_for_DG(fieldname, value)[source]

Returns the weak Dirichlet terms for a discontinuous Galerkin (DG) formulation. When using a DirichletBC with prefer_weak_for_DG, this method is called. If it returns not None, the DirichletBC is not enforced strongly, but on the basis of the given interface residuals.

Parameters:

fieldname (str) – Name of the field for which the weak Dirichlet terms are to be returned.
value (Union[Expression, int, float]) – The desired Dirichlet condition.

Return type:

Union[Expression, int, float, None]

class pyoomph.generic.GenericProblemHooks[source]

Bases: object

A class that can be attached to a problem to call additional functions after e.g. newton solves, etc.

class pyoomph.generic.GlobalLagrangeMultiplier(*args, only_for_stationary_solve=False, set_zero_on_normal_mode_eigensolve=True, **kwargs)[source]

Bases: ODEEquations

This class allows to define a global Lagrange multiplier that are used to enforce global constraints. It is just a normal ODEEquations but with some additional features, i.e. it can be e.g. deactivated on transient solves.

class pyoomph.generic.InterfaceEquations(*args: Any, **kwargs: Any)[source]

Bases: Equations

Same as normal Equations but with some extra functions for equations defined on interfaces.

get_opposite_parent_equations(of_type=None)[source]

Returns the Equations in the parent bulk domain at the opposite side of this interface. When setting the attribute required_opposite_parent_type, of_type can be omitted to get the expected parent equations.

This method is useful to e.g. get the mass density from a Navier-Stokes equation in the opposite bulk domain, e.g. for mass transfer processes at the interface.

Parameters:: of_type (Optional[Type[Equations]]) – The type of the equations to be returned. If not provided, the required_opposite_parent_type has to be set.

get_parent_equations(of_type=None)[source]

Returns the Equations in the parent bulk domain of a given type. When setting the attribute required_parent_type, of_type can be omitted to get the expected parent equations.

This method is useful to e.g. get the mass density from a Navier-Stokes equation in the bulk domain, e.g. for mass transfer processes at the interface.

Parameters:: of_type (Optional[Type[Equations]]) – The type of the equations to be returned. If not provided, the required_parent_type has to be set.

pin_redundant_lagrange_multipliers(mesh, lagr, depvars, opposite_interface=[])[source]

Allows to pin redundant (overconstraining) Lagrange multipliers. A field of Lagrange multipliers usually enforces some constraint depending on depvars (and poentially degrees at the opposite_interface). If all these degrres are pinned, the Lagrange multiplier lagr is pinned and set to zero as well.

Parameters:

mesh (InterfaceMesh) – The current mesh must be passed
lagr (str) – Name of the Lagrange multiplier field to be automatically pinned if necessary
depvars (Union[str, List[str], Tuple[str, ...]]) – Single or multiple variables that occur in the constraint.
opposite_interface (Union[str, List[str], Tuple[str, ...]]) – Optional dependencies on the opposite side of the interface.

required_opposite_parent_type: Optional[Type[Equations]] = None: If set to a particular Equations class, pyoomph will check whether we have indeed these equations at the opposite bulk side of this interface

required_parent_type: Optional[Type[Equations]] = None: If set to a particular Equations class, pyoomph will check whether we have indeed these equations in the bulk

class pyoomph.generic.ODEEquations(*args: Any, **kwargs: Any)[source]

Bases: BaseEquations

Class representing a set of ordinary differential equations (ODEs). Add unknowns by overriding the define_fields() method and residuals by overriding the define_residuals() method.

See BaseEquations for further methods.

define_ode_variable(*names, scale=None, testscale=None)[source]

Defines the ODE variables.

Parameters:

*names (str) – Variable length argument representing the names of the ODE variable(s).
scale (Union[Expression, int, float, None]) – Optional scaling factor for the ODE variable(s).
testscale (Union[Expression, int, float, None]) – Optional scaling factor for the test functions associated with the ODE variable(s).

Return type:

None

expand_additional_field(name, dimensional, expression, in_domain, no_jacobian, no_hessian, where)[source]

Expands additional fields for the ODEs.

Parameters:

name (str) – The name of the additional field.
dimensional (bool) – A boolean indicating if the field is dimensional.
expression (Expression) – The expression defining the additional field.
in_domain (FiniteElementCode) – The finite element code representing the domain.
no_jacobian (bool) – A boolean indicating if the Jacobian should not be computed.
no_hessian (bool) – A boolean indicating if the Hessian should not be computed.
where (str) – The location where the additional field is expanded.

Returns:

The expanded additional field.

Return type:

Expression

get_coordinate_system()[source]

Returns the base coordinate system.

Returns:: The base coordinate system.
Return type:: BaseCoordinateSystem

get_mesh()[source]

Returns the ODE storage mesh.

Returns:: The ODE storage mesh.
Return type:: ODEStorageMesh

get_parent_domain()[source]

Returns the parent domain.

Returns:: The parent domain.
Return type:: None

on_apply_boundary_conditions(mesh)[source]

Applies boundary conditions to the ODEs.

Parameters:: mesh (Union[InterfaceMesh, MeshFromTemplate1d, MeshFromTemplate2d, MeshFromTemplate3d, ODEStorageMesh]) – The mesh to which the boundary conditions are applied.

pin(**kwargs)[source]

Pins the specified degrees of freedom (DOFs).

Parameters:: **kwargs (Union[bool, float]) – Keyword arguments representing the DOFs to be pinned and their values. If assigning True, we just fix the current value, otherwise, the provided value is set.

unpin(*args)[source]

Unpins the specified degrees of freedom (DOFs).

Parameters:: *args (str) – Variable length argument representing the DOFs to be unpinned.

class pyoomph.generic.Problem[source]

Bases: Problem

A class representing a problem in the pyoomph library.

This class provides methods and attributes for defining and solving a problem. Usually, in the __init__() method, you define any parameters with default settings. The problem itself is defined by the define_problem() method, where you define the equations and the mesh(es). After creation of an instance, you can solve the problem by calling the run() method (transient solves) or the solve() method (stationary solve). Outputs (potentially with plots) can be generated by calling the output() method.

additional_equations

Additional equations for the problem.

Type:: Union[Literal[0], EquationTree]

always_take_one_newton_step

Flag indicating whether to always take one Newton step.

Type:: bool

continuation_data_in_states

Flag indicating whether to store continuation data in the states.

Type:: bool

default_1d_file_extension

Default file extension for 1D files.

Type:: Union[Literal[“txt”, “mat”], List[Literal[“txt”, “mat”]]]

default_ccode_expression_mode

Default C code expression mode.

Type:: str

default_spatial_integration_order

Default spatial integration order.

Type:: Union[int, None]

default_timestepping_scheme

Default timestepping scheme.

Type:: Literal[“BDF2”, “BDF1”, “Newmark2”]

eigen_data_in_states

Flag indicating whether to store eigen data in the states.

Type:: Union[int, bool]

eigenvector_position_scale

Scaling factor for eigenvector positions.

Type:: float

extra_compiler_flags

Extra compiler flags for the problem.

Type:: List[str]

ignore_command_line

Flag indicating whether to ignore command line arguments.

Type:: bool

latex_printer

LaTeX printer for the problem.

Type:: Optional[LaTeXPrinter]

max_residuals

Maximum residuals for the problem.

Type:: float

plot_in_dedicated_process

Flag indicating whether to plot in a dedicated process.

Type:: bool

remove_macro_elements_after_initial_adaption

Flag indicating whether to remove macro elements after initial adaption.

Type:: Union[bool, Literal[“auto”]]

scaling

Dictionary of scaling factors.

Type:: Dict[str, Union[str, ExpressionOrNum]]

states_compression_level

Compression level for the states.

Type:: Union[int, None]

timestepper

Timestepper for the problem.

Type:: MultiTimeStepper

write_states

Flag indicating whether to write states.

Type:: bool

abort_current_run()[source]: If called within a run(…) statement, e.g. from some action_{after/before}_* methods, the run will abort

activate_bifurcation_tracking(parameter, bifurcation_type=None, blocksolve=False, eigenvector=None, omega=None, azimuthal_mode=None, cartesian_wavenumber_k=None, eigenvalue_for_branch_tracking=None)[source]

Activates bifurcation tracking for the specified parameter and bifurcation type. Subsequent calls of solve(…) and arclength_continuation(…) will then track the bifurcation.

Parameters:

parameter (Union[str, GiNaC_GlobalParam, None]) – The parameter to change in order to find the bifurcation. If None, we track the current eigenbranch, i.e. Re(lambda) will be found and is not necessarily 0.
bifurcation_type (Optional[Literal["hopf", "fold", "pitchfork", "azimuthal"]]) – The type of bifurcation to track. Defaults to None, i.e. auto-detect.
blocksolve (bool) – Flag indicating whether to use block solve. Defaults to False. Should be kept False.
eigenvector (Optional[Union[NPFloatArray, NPComplexArray, int]]) – The eigenvector to use for tracking. Defaults to None, which means the eigenvector corresponding to the eigenvalue with largest real part. Can be either an index or a custom vector.
omega (Optional[float]) – The omega value for Hopf bifurcation tracking. Defaults to None, then it will be Im(lambda).
azimuthal_mode (Optional[int]) – The azimuthal mode for azimuthal bifurcation tracking. Defaults to None.

activate_eigenbranch_tracking(branch_type=None, eigenvector=None, eigenvalue=None)[source]

Activates eigenbranch tracking for the specified eigenbranch type. Subsequent calls of solve(…) and arclength_continuation(…) will then track the eigenbranch. This is similar to bifurcation tracking, but it does not adjust a parameter to find a bifurcation, i.e. where Re(lambda)=0. Instead, it starts with a eigenvalue/eigenvector pair. Once activated, you can follow the eigenbranch by calling arclength_continuation(…). At each step, the eigenvalue/eigenvector pair will be updated and is available via get_last_eigenvalues()[0] and get_last_eigenvectors()[0].

Parameters:

branch_type (Optional[Literal["real", "complex", "normal_mode"]]) – The type of eigenbranch to track. Defaults to None, i.e. auto-detect.
eigenvector (Optional[int]) – The previously calculated eigenvector index to use for tracking. Defaults to None, i.e. the eigenvector at index zero.

activate_periodic_orbit_handler(T, history_dofs=[], mode='collocation', order=2, GL_order=-1, T_constraint='phase')[source]

Activates periodic orbit tracking based on history dofs. Use set_current_dofs() to set the first time point of the orbit guess. The other time points must be shipped with the history_dofs argument.

Parameters:

T (Union[Expression, int, float]) – The guessed period of the orbit
history_dofs – The history dofs to use for the orbit tracking. Must be non-empty.
mode (Literal['collocation', 'floquet', 'bspline', 'central', 'BDF2']) – The mode of the time discretization.
order (int) – The order of the time discretization.
GL_order (int) – The Gauss-Legendre order for some time discretization modes. Defaults to -1, meaning a suitable integration order is chosen automatically based on the interpolation order.
T_constraint (Literal['plane', 'phase']) – The constraint for the period. Defaults to “phase”.

Returns:

The resulting periodic orbit. Note that it still must be solved, i.e. it is only the provided guess at this stage.

Return type:

PeriodicOrbit

add_equations(eqs)[source]

Add equations to the system. Should be called within the define_problem() method.

Parameters:: eqs (EquationTree) – The equations (restricted to a domain) to add to the problem.
Return type:: None

add_global_dof(name, equation_contribution=0, *, scaling=None, testscaling=None, domain='globals', only_for_stationary_solve=False, initial_condition=None, set_zero_on_normal_mode_eigensolve=True)[source]

Add a global degree of freedom, e.g. a global Lagrange multiplier to the problem.

Parameters:

name (str) – The name of the degree of freedom.
equation_contribution (ExpressionOrNum) – The global contribution of the degree of freedom to its equation. Defaults to 0.
scaling (ExpressionNumOrNone, optional) – The scaling factor for the degree of freedom. Defaults to None.
testscaling (ExpressionNumOrNone, optional) – The scaling factor for the test function. Defaults to None.
domain (str, optional) – The domain to which the degree of freedom belongs. Defaults to “globals”.
only_for_stationary_solve (bool, optional) – Whether the degree of freedom is only used for stationary solves, if set, it will be 0 and pinned during transient solves. Defaults to False.
initial_condition (ExpressionNumOrNone, optional) – The initial condition for the degree of freedom. Defaults to None.
set_zero_on_normal_mode_eigensolve (bool) – Deactivate this dof for normal mode eigensolves. Defaults to True.

Returns:

A tuple containing the variable and test function associated with the degree of freedom.

Return type:

tuple

add_mesh(mesh)[source]

Adds a mesh to the problem. Based on the domain and boundary names of the mesh, equations can be added by using the same domain and boundary names.

Parameters:: mesh (TypeVar(_TypeVarMeshTemplate, bound= MeshTemplate)) – Any mesh instance to be added (1d, 2d, 3d, etc.)
Return type:: TypeVar(_TypeVarMeshTemplate, bound= MeshTemplate)
Returns:: Returns itself for chaining

add_mesh_template(mesh)[source]

Same as self+=mesh or self.add_mesh(mesh). Will be deprecated soon.

Return type:: TypeVar(_TypeVarMeshTemplate, bound= MeshTemplate)

arclength_continuation(parameter, step, *, spatial_adapt=0, max_ds=None, max_newton_iterations=None, newton_relaxation_factor=None, newton_solver_tolerance=None, min_ds=None, dof_direction=None, globally_convergent_newton=False)[source]

Perform arclength continuation on the basis of a given parameter.

Parameters:

parameter (Union[str, _pyoomph.GiNaC_GlobalParam]) – The parameter to perform arclength continuation on.
step (float) – The step for the continuation.
spatial_adapt (int, optional) – The level of spatial adaptation. Defaults to 0.
max_ds (float, optional) – The maximum step size. Defaults to None.
max_newton_iterations (int, optional) – The maximum number of Newton iterations. Defaults to None.
newton_relaxation_factor (float, optional) – The relaxation factor for the Newton solver. Defaults to None.
newton_solver_tolerance (float, optional) – The tolerance for the Newton solver. Defaults to None.
min_ds (float, optional) – The minimum step size. Defaults to None.
dof_direction (List[float], optional) – The direction of degrees of freedom. Defaults to None.
globally_convergent_newton (bool, optional) – Whether to use globally convergent Newton solver. Defaults to False.

Returns:

The new step size for the continuation.

Return type:

float

check_mesh_integrity: Union[bool, Literal['initially']]: Checks whether the elements in the meshes are nicely oriented (facing) so that refinement works as it should. Can be only done once initially or at each refinement step

continue_from_outdir(old_out_dir, statenumber=-1, ignore_outstep=True)[source]

Loads a previous state from another output directory. Make sure the scripts are in all specifications of equations, meshes, parameters, settings etc.

Parameters:

old_out_dir (str) – Old output directory
statenumber (int) – Which state file to load (default: -1, i.e. the last one)
ignore_outstep (bool) – Do not load the outstep (default: True)

create_eigendynamics_animation(outdir, plotter, eigenvector=0, init_amplitude=None, max_amplitude=None, numperiods=1, numouts=25, phi0=0)[source]

Creates an animation of the eigenfunction dynamics. The eigenfunction is animated by varying the time and the amplitude of the eigenfunction, which is added to the degrees of freedom at each time. All images are saved in the specified output directory (relative to the output directory of the problem). The plotter is used to create the images. Azimuthal instabilities will automatically mirror the eigenfunction to the left in the appropriate way.

Parameters:

outdir (str) – Output directory for the animation images relative to the output directory of the problem.
plotter (MatplotlibPlotter) – Plotter class to use for the animation.
eigenvector (int) – Optional index of the eigenfunction to animate. Defaults to 0.
init_amplitude (Optional[float]) – Initial amplitude of the eigenperturbation. If this and max_amplitude is not provided, the amplitude is set to 1. Defaults to None.
max_amplitude (Optional[float]) – Maximum amplitude of the eigenperturbation. If this is provided, the amplitude is set to this value at the beginning and decreases over time (eigenvalue has negative real part) or will reach this amplitude at the end of the considered time (eigenvalue has positive real part). Defaults to None.
numperiods (float) – Number of periods to animate. For purely real eigenvalues, the characteristic time is given by the real part. Defaults to 1.
numouts (int) – Number of output steps. Defaults to 25.
phi0 (float) – Initial phase. Defaults to 0.

create_text_file_output(filename, header=None, relative_to_output_dir=True)[source]

Creates a NumericalTextOutputFile. By default, in the output directory.

Parameters:

filename (str) – File name
header (Optional[List[str]]) – Header of the file Defaults to None.
relative_to_output_dir (bool) – If True, the file is created in the output directory. Defaults to True.

Return type:

NumericalTextOutputFile

deactivate_bifurcation_tracking()[source]: Deactivate bifurcation tracking. Afterwards, the problem can be solved as usual.

define_global_parameter(**params)[source]

Define one or more global parameters for the problem.

Parameters:: **params (float) – A dictionary of parameter names and their corresponding initial values.
Returns:: The defined global parameter(s), which can be used in expressions.
Return type:: Union[_pyoomph.GiNaC_GlobalParam, Tuple[_pyoomph.GiNaC_GlobalParam, …]]

define_named_var(**kwargs)[source]

Named vars are global expressions that are bound to a name.

You can e.g. use define_named_var(temperature=20*celsius) to define a temperature variable at problem level. When an equation tried to expand var(“temperature”) and no field “temperature” is defined on the current domain or its parents, the variable will be expanded by the global variable.

define_problem()[source]

Define the problem by creating the mesh(es) and other necessary components.

This method should be overridden by subclasses to define the specific problem.

deflated_continuation(deflation_alpha=0.1, deflation_p=2, perturbation_amplitude=0.5, max_newton_iterations=None, newton_relaxation_factor=None, use_eigenperturbation=False, skip_initial_solution=False, num_random_tries=1, max_branches=None, branch_continue_iterations=10, **param_range)[source]

Scan over a parameter range and try to find multiple solutions for each parameter step by deflation This is an implemetation according to: The computation of disconnected bifurcation diagrams by Patrick E. Farrell, Casper H. L. Beentjes, Ásgeir Birkisson

https://arxiv.org/pdf/1603.00809.pdf

Parameters:

deflation_alpha (float) – Shift of the deflation operator. Defaults to 0.1.
deflation_p (int) – Order of the deflation. Defaults to 2.
perturbation_amplitude (float) – Perturbation amplitude to move away from the previous solution. Defaults to 0.5.
max_newton_iterations (Optional[int]) – Optional override of the number of Newton iterations during deflated search for additional solutions. Defaults to None.
newton_relaxation_factor (Optional[float]) – Optional override of the Newton relaxation factor during deflated search for additional solutions. Defaults to None.
use_eigenperturbation (bool) – Whether to use eigen perturbation for the next solution during deflation. Defaults to False.
num_random_tries (int) – Number of random tries for finding solutions during deflation. Defaults to 1.
max_branches (Optional[int]) – Maximum number of branches to find. Defaults to None.
branch_continue_iterations (int) – Number of iterations for continuing branches. Defaults to 10.

Yields:

A tuple of branch index (from 0 to …), the current parameter value and the current degrees of freedom (dofs) for the solution.

deflated_solve_by_eigenperturbation(eigenindex=0, keep_deflation_active=False, perturbation_factor=1, deflation_alpha=0.1, deflation_power=2, *, max_newton_iterations=None, newton_relaxation_factor=None, newton_solver_tolerance=None, globally_convergent_newton=False)[source]

Tries to find another stationary solution by deflation. The procedure is implemented according to ‘Deflation techniques for finding distinct solutions of nonlinear partial differential equations’ by Patrick E. Farrell, Ásgeir Birkisson & Simon W. Funke, https://arxiv.org/pdf/1410.5620.pdf .

Args:
deflation_alpha (float, optional): Shift of the deflation operator. Defaults to 0.1. deflation_p (int, optional): Order of the deflation. Defaults to 2. perturbation_amplitude (float, optional): Perturbation amplitude to move away from the previous solution. Defaults to 1. max_newton_iterations (Optional[int], optional): Optional override of the number of Newton iterations to try. Defaults to None. newton_relaxation_factor (Optional[float], optional): Optional override of the Newton relaxation factor. Defaults to None.

do_call_remeshing_when_necessary: bool: Flag indicating whether to call remeshing when necessary. Can be set to False to disable remeshing, e.g. when tracking bifurcations, it is better to check it manually invoking the remeshing method remesh_handler_during_continuation() after each solve.

dof_strings_to_global_equations(string_dof_set)[source]

Takes strings like "domain/velocity_x" and returns a set of global equations

Parameters:: string_dof_set (Union[str, Set[str], List[str]]) – Degrees of freedom you want to resolve to equation numbers
Returns:: Global equation set
Return type:: Set[int]

find_bifurcation_via_eigenvalues(parameter, initstep, shift=0, neigen=6, spatial_adapt=0, epsilon=1e-08, reset_arclength=False, max_ds=None, stay_stable_file=None, before_eigensolving=None, do_solve=True, azimuthal_m=None, normal_mode_k=None, eigenindex=0)[source]

Approximates a bifurcation point by bisecting on the basis of the eigenvalues. Must be called as a generator, e.g.

for parameter, eigenvalue in find_bifurcation_via_eigenvalues(...):
    print("Currently at ",parameter,eigenvalue)

Parameters:

parameter (Union[str,_pyoomph.GiNaC_GlobalParam]) – The parameter to vary to find a bifurcation. It can be either a string representing the name of a global parameter or a global parameter directly.
initstep (float) – The initial step size for the bisection.
shift (Union[None,float,complex]) – The shift value for the eigenvalue problem. It can be a float, complex number, or None.
neigen (int) – The number of eigenvalues to compute.
spatial_adapt (int) – The spatial adaptation level.
epsilon (float) – The tolerance for determining the real part of an eigenvalue to be close to zero.
reset_arclength (bool) – Whether to reset the arc length parameters after each step.
max_ds (Optional[Union[float,Callable[[float],float]]]) – The maximum step size for the continuation. It can be a float, a callable function that takes the current parameter value and returns a float, or None.
stay_stable_file (str, optional) – The file path to save the state when the solution is stable. If it is unstable, the state is reloaded.
before_eigensolving (Optional[Callable[[float],None]]) – A callable function to be called before solving the eigenvalue problem.
do_solve (bool) – Whether to solve the problem before continuation. If the solution does not depend on the parameter, it can be set to False.
azimuthal_m (Union[int, List[int], None]) – The azimuthal mode number if you want to find azimuthal perturbations.

Yields:

param(float) – The current parameter value. eigenvalue (complex): The eigenvalue corresponding to the specified eigenindex.

Raises:

RuntimeError – If the eigenindex is greater than or equal to neigen.
RuntimeError – If the initial solution is already unstable.

get_cached_mesh_data(msh, nondimensional=False, tesselate_tri=True, eigenvector=None, eigenmode='abs', history_index=0, with_halos=False, operator=None, discontinuous=False, add_eigen_to_mesh_positions=True)[source]

Return the current data (i.e. values) of a mesh. These are cached in case they are required multiple times, e.g. for plotting and output. The cache is invalidated whenever we solve the problem or set some initial condition.

Parameters:

msh (Union[str,AnySpatialMesh]) – Mesh object or mesh name
nondimensional (bool, optional) – Getting nondimensional values instead of dimensional ones. Defaults to False.
tesselate_tri (bool, optional) – Split quad elements into tris. Helpful e.g. for plotting via matplotlib, which requires triangular meshes. Defaults to True.
eigenvector (Optional[Union[int,Sequence[int]]], optional) – If not None, we can obtain the values of the eigenfunction with the given index. Defaults to None.
eigenmode (MeshDataEigenModes, optional) – Since eigenfunctions are in general complex, we must select the desired projection to real numbers here in case of eigenvector!=None. Defaults to “abs”.
history_index (int, optional) – Set to 1 or 2 to access the previous time steps. Defaults to 0.
with_halos (bool, optional) – Include halos to the output. Defaults to False.
operator (Optional[MeshDataCacheOperatorBase], optional) – Apply an operator on the cache, e.g. to add eigenvectors or extrude it in 3d. Defaults to None.

Returns:

The combined information of the mesh cache

Return type:

MeshDataCacheEntry

get_coordinate_system()[source]

Get the coordinate system set at problem level.

Returns:: The coordinate system at problem level.
Return type:: BaseCoordinateSystem

get_current_time(dimensional=True, as_float=False)[source]

Get the current time of the problem.

Parameters:

dimensional (bool, optional) – Flag indicating whether to return the dimensional time. Defaults to True.
as_float (bool, optional) – Flag indicating whether to return the time as a float. Defaults to False.

Returns:

The current time of the problem.

Return type:

ExpressionOrNum

Raises:

ValueError – If the problem is not initialized.

get_dof_description()[source]

Returns two arrays containing the description of the degrees of freedom. The first is a list of dof-type indices, where the i-th entry is the type of the i-th degree of freedom. To resolve what each dof-type index means, the second array contains the type names of the degrees of freedom.

For a simple Poisson equation for a field u on a line domain name "domain", the first returned array will contain numbers between 0 and 2 (dof-type indices), and the second array will contain the type names ["domain/u", "domain/left/u", "domain/right/u"]. Note how the boundaries get their own dof-type indices.

Returns:: A pair of arrays containing the dof-type indices and the type names to classify the degrees of freedom.

get_eigen_solver()[source]

Get the eigenproblem solver instance.

Returns:: The currently used eigensolver
Return type:: GenericEigenSolver

get_equations(path, error_if_not_found=True)[source]

Return the equations added at the specified path.

Parameters:

path (str) – Path to the domain
error_if_not_found (bool, optional) – Raise an error if the domain path is not valid. Defaults to True.

Raises:

RuntimeError – In case you specify a path where not equations are defined and error_if_not_found==True, you will get this error

Returns:

The equations at the desired domain path. If error_if_not_found==False, it is None if there are no equations at the desired path specified.

Return type:

Optional[BaseEquations]

get_floquet_multipliers(n=None, valid_threshold=10000, shift=None, ignore_periodic_unity=False, quiet=True)[source]

TODO; Add documentation

Return type:: ndarray[tuple[Any, ...], dtype[complex128]]

get_la_solver()[source]

Get the linear solver instance.

Returns:: The currently used linear solver
Return type:: GenericLinearSystemSolver

get_last_eigenmodes_k()[source]

Get the cartesian normal mode numbers for the last computed eigenvalues.

Returns:: Array containing the cartesian normal mode numbers corresponding to the eigenvalues.
Return type:: Optional[NPFloatArray]

get_last_eigenmodes_m()[source]

Get the azimuthal mode numbers for the last computed eigenvalues.

Returns:: Array containing the azimuthal mode numbers corresponding to the eigenvalues.
Return type:: Optional[NPIntArray]

get_last_eigenvalues(dimensional=False)[source]

Returns the last computed eigenvalues.

Returns:: Eigenvalues as array.
Return type:: NPComplexArray

get_last_eigenvectors()[source]

Return the last computed eigenvector.

Returns:: Eigenvectors as 2d array.
Return type:: NPComplexArray

get_mesh(name, return_None_if_not_found=False)[source]

Get the mesh at the desired domain path. Invokes initialization if the problem is not initialised!

Parameters:

name (str) – Domain path of the mesh
return_None_if_not_found (bool, optional) – If True, None will be returned if the given domain path is invalid. If False, an error will be raised in that case. Defaults to False.

Raises:

RuntimeError – Raised if there is no mesh at the given domain path in case of return_None_if_not_found==False (default). Same happens if an ODE domain is tried to be accessed like this. Use get_ode() for this.

Returns:

The mesh at the domain path. None can only be returned if return_None_if_not_found==True and the domain path is invalid.

Return type:

Optional[AnySpatialMesh]

get_ode(name)[source]

Return the ODE object at the given domain path. Invokes initialization if the problem is not initialised!

Parameters:: name (str) – Domain path of the ODE
Raises:: RuntimeError – If the given domain path is invalid or a spatial mesh is defined at this domain, this error will occur.
Returns:: The ODE object at the given domain path
Return type:: ODEStorageMesh

get_output_directory(relative_path=None)[source]

Return the output directory of the problem. Set it with set_output_directory(). Otherwise, it will default to the name of the invoked script minus the extension .py. Optionally, you can add a relative path to assemble e.g. a file name within the output directory.

Parameters:: relative_path (Optional[str], optional) – If set, we join this relative additional path to the output directory. Defaults to None.
Returns:: The output directory of the problem (potentially joined with the additionally passed relative_path)
Return type:: str

get_scaling(s, none_if_not_set=False)[source]

Get the scaling factor for the problem variables for nondimensionalization.

Parameters:

s (str) – Name of the scale to get.
none_if_not_set (bool) – Returns None if this scaling is not set. Otherwise, the default scale 1 is returned. Defaults to False.

Return type:

Union[Expression, int, float, None]

Returns:

Scaling set by set_scaling() or None if none_if_not_set==True and the scale is not set.

gitignore_output: bool: Add a .gitignore with content “*” to output folders

go_to_param(*, reset_pars=True, startstep=None, call_after_step=None, final_adaptive_solve=False, max_newton_iterations=None, epsilon=1e-06, max_step=None, **kwargs)[source]

Perform arclength continuation in a parameter until we reach the desired value.

Parameters:

reset_pars (bool, optional) – Whether to reset arc length parameters. Defaults to True.
startstep (float, optional) – The initial step size for the parameter continuation. Defaults to None.
call_after_step (Callable[[float],None], optional) – A function to call after each step. If it returns “stop”, we stop any further continuation. Defaults to None.
final_adaptive_solve (Union[bool,int], optional) – Whether to perform a final adaptive solve. Defaults to False.
max_newton_iterations (int, optional) – The maximum number of Newton iterations. Defaults to None.
epsilon (float) – The tolerance for considering as converged to the parameter
max_step (Optional[float]) – The maximum step size for the continuation. Defaults to None.
**kwargs (float) – The parameter name and desired value.

Raises:

RuntimeError – If more than one parameter is provided.
RuntimeError – If the specified parameter is not part of the problem.

Return type:

None

Returns:

None

guess_nearest_bifurcation_type(eigenvector=0)[source]: Guesses the nearest bifurcation type based on the last computed eigenvalues. This is only possible after calling solve_eigenproblem(…). It cannot guess e.g. pitchfork or transcritical bifurcations, only “hopf” or “fold” - or “azimuthal” if the last eigenvalues correspond to azimuthal modes m!=0. :returns: Guessed bifurcation type :rtype: str

initial_adaption_steps: Optional[int]: Spatial adaption steps for the initial condition. If set to None, we refine initially up to max_refinement_level.

initialise()[source]

Initializes the problem by performing the necessary setup and initialization steps. If not done before, this method is automatically called by several methods, e.g. solve(), run() or output(). After initialization, you cannot change the problem anymore, except for global parameter values.

Raises:: RuntimeError – If the problem is already initialised or if a function that calls initialize is called during initialization.

invalidate_cached_mesh_data(only_eigens=False)[source]

Mesh data is cached for potentially multiple usage (e.g. plotting and output to file). Whenever we change anything (e.g. changing values), we must hence invalidate the cache.

Parameters:: only_eigens (bool, optional) – Only flush the cache of eigenfunctions, not of the current state. Defaults to False.

is_initialised()[source]

Returns whether the problem has been initialised or not.

Returns:: True if already initialised, False otherwise
Return type:: bool

is_normal_mode_stability_set_up()[source]

Returns True when setup_for_stability_analysis() has been called with azimuthal_stability=True or additional_cartesian_mode=True. Can be used to e.g. set additional BCs for velocity_phi or similar.

Return type:: Union[Literal['azimuthal', 'cartesian'], Literal[False]]

is_stable_solution()[source]

Shortcut to check whether we have a stable solution. This is only possible after calling solve_eigenproblem(…).

Return type:: bool

iterate_over_multiple_solutions_by_deflation(deflation_alpha=0.1, deflation_p=2, perturbation_amplitude=0.5, max_newton_iterations=None, newton_relaxation_factor=None, use_eigenperturbation=False, skip_initial_solution=False, num_random_tries=1, keep_deflation_operator_active=False)[source]

Tries to find multiple stationary solutions by deflation. The procedure is implemented according to ‘Deflation techniques for finding distinct solutions of nonlinear partial differential equations’ by Patrick E. Farrell, Ásgeir Birkisson & Simon W. Funke, https://arxiv.org/pdf/1410.5620.pdf . :rtype: Generator[ndarray[tuple[Any, ...], dtype[float64]], None, None]

Args:
deflation_alpha (float, optional): Shift of the deflation operator. Defaults to 0.1. deflation_p (int, optional): Order of the deflation. Defaults to 2. perturbation_amplitude (float, optional): Perturbation amplitude to move away from the previous solution. Defaults to 0.5. max_newton_iterations (Optional[int], optional): Optional override of the number of Newton iterations to try. Defaults to None. newton_relaxation_factor (Optional[float], optional): Optional override of the Newton relaxation factor. Defaults to None.

Yields:
The found solutions as lists of degrees of freedom

logfile_name: Optional[str]: Name of the logfile (or None for no logfile), relative to the output directory

max_permitted_error: float: Maximum error of all meshes for spatial adaptivity. If the error is above this threshold, we must refine locally.

max_refinement_level: int: Maximum number of refinements of all meshes.

min_permitted_error: float: Minimum error of all meshes for spatial adaptivity. If the error is below this threshold, we may unrefine locally.

min_refinement_level: int: Minimum refinement level of all meshes.

output(stage='', quiet=None)[source]

Invoke an output of the current solution at the current time by calling all Output objects.

Parameters:

stage (str) – The stage of the output, at the moment, only “” is meaninfull.
quiet (bool, optional) – Flag to control the verbosity of the output.

Return type:

None

Returns:

None

output_at_increased_time(dt=None)[source]

Increases the current time by the specified time step (dt, default scale_factor(“temporal”)) and calls the output method. Useful for Paraview PVD output of multiple stationary solutions, which otherwise overlays multiple outputs at the same time step.

Parameters:: dt (Optional[ExpressionOrNum]) – The time step to increase the current time by. If not provided, the scaling factor for temporal is used.
Return type:: None
Returns:: None

perturb_dofs(dofpert)[source]

Perturbs all degrees of freedom by a given perturbation array (must have the length of ndof())

Parameters:: dofpert (ndarray[tuple[Any, ...], dtype[float64]]) – Perturbation array to be added to the degrees of freedom (nondimensional)

plotter: Union[List[MatplotlibPlotter], MatplotlibPlotter, None]: Set this to a (list of ) plotter(s) to automatically plot on output() calls. If set to None, no plotting will be done.

precice_config_file: str: Must be set to the config file when using preCICE

precice_initialise()[source]: Initializes the preCICE adapter for the problem. You must set precice_participant and precice_config_file in the Problem class.

precice_participant: str: Must be set to the participant name when using preCICE. Default is an empty string, if you do not use preCICE.

precice_run(maxstep=None, temporal_error=None, output_initially=True, fast_dof_backup=False)[source]

Runs a simulation with the precice adapter. To that end, you must set precice_participant and precice_config_file in the Problem class. There is less control compared to the normal py:meth:pyoomph.generic.problem.Problem.run (i.e. without preCICE), but a lot of settings can be adjusted in the preCICE configuration file.

Parameters:

maxstep (Optional[float]) – Maximum nondimensional time step. Defaults to None.
temporal_error (Optional[float]) – Use temporal adaptivity with this given error factor. Defaults to None.
output_initially (bool) – Outputs before the simulation starts. Defaults to True.
fast_dof_backup (bool) – If True, only the DoFs will be backed up, nothing else. Defaults to False.

process_eigenvectors(eigenvectors)[source]

Here, you can optionally scale or negate the eigenvectors, e.g. normalize them, by overriding this method

Parameters:: eigenvectors (NPFloatArray) – 2d array of eigenvectors
Returns:: Processed array of eigenvectors
Return type:: NPFloatArray

redefine_problem(code_dir, interpolator=<class 'pyoomph.meshes.interpolator.InternalInterpolator'>, num_adapt=None)[source]

Redefines the problem by recompiling equations. This can in principle be used if problem parameters have changed, but it is not recommended to change the problem structure. If possible, it is advised to use the global parameter system to change any parameters.

Parameters:

code_dir (str) – Subdirectory in the output directory where the C++ code of the redefined problem will be written.
interpolator (Type[BaseMeshToMeshInterpolator]) – Mesh interpolator to map the fields of the old meshes to the new ones.
num_adapt (Optional[int]) – Number of adaption steps after redefining the problem. If None, the number of adaption steps is determined by the max_refinement_level attribute.

Raises:

RuntimeError – If the problem contains no equations after the redefinition

refine_eigenfunction(numadapt=1, eigenindex=0, resolve_base_state=True, resolve_neigen=1, use_startvector=False)[source]

After calculating an eigenproblem, you can adapt the mesh according ot the eigenfunction of a specific eigenvalue. This can be useful to refine the mesh in regions where the eigenfunction has a high gradient. It requires SpatialErrorEsitmators to be added to the problem and an adaptive mesh.

Parameters:

numadapt (int) – The number of adaptations to perform. Defaults to 1.
eigenindex (int) – The index of the eigenvalue to refine the mesh for. Defaults to 0.
resolve_base_state (bool) – If True, the base state is resolved after each adaptation. Defaults to True.
resolve_neigen (int) – The number of eigenvalues to resolve after each adaptation. Defaults to 1.

Returns:

The eigenvalue and eigenvector of the adapted eigenproblem.

Return type:

Tuple[float, NPFloatArray]

remesh_handler_during_continuation(force=False, resolve=True, resolve_before_eigen=False, reactivate_biftrack_neigen=4, reactivate_biftrack_shift=0, resolve_max_newton_steps=None, num_adapt=None, resolve_globally_convergent_newton=False)[source]

Handle remeshing during continuation. We might have to calculate e.g. a new eigenvector when doing bifurcation tracking. In that case, set Problem.do_remeshing_when_necessary to False to prevent any automatic remeshing.

Parameters:

force (bool, optional) – Force remeshing even if not necessary. Defaults to False.
resolve (bool, optional) – Resolve the problem after remeshing. Defaults to True.
resolve_before_eigen (bool, optional) – Resolve the problem before solving the eigenproblem. Defaults to False.
reactivate_biftrack_neigen (int, optional) – Number of eigenvalues to reactivate bifurcation tracking. Defaults to 4.
reactivate_biftrack_shift (float, optional) – Shift for the eigenvalues to reactivate bifurcation tracking. Defaults to 0.
resolve_max_newton_steps (int, optional) – Maximum number of Newton steps to resolve the problem.
resolve_globally_convergent_newton (bool) – Use a globally convergent Newton solver. Defaults to False.

Returns:

True if remeshing was performed, False otherwise.

Return type:

bool

remesh_if_necessary()[source]

Invokes remeshing if one RemeshWhen object indicates that.

Returns:: True if remeshing was performed, False otherwise.
Return type:: bool

remeshing_necessary()[source]

Checks whether any RemeshWhen object indicates that remeshing should be done.

Returns:: True if remeshing would be required, False otherwise.
Return type:: bool

rotate_eigenvectors(eigenvectors, dofs_to_real, normalize_dofs=False, normalize_amplitude=1, normalize_max=True)[source]

Should be called within the method process_eigenvectors() to rotate the eigenvectors to e.g. a common phase. This is optional, but avoids phase jumps in the eigenvectors when following an eigenbranch.

Parameters:

eigenvectors – Eigenvectors to rotate, usually the ones passed in the automatically method process_eigenvectors().
dofs_to_real (Union[str, List[str], Set[str]]) – Which degrees of freedom to consider to find the phase. Can be a single string, a list of strings or a set of strings.
normalize_dofs (bool) – Normalizes the eigenvector with respect to the selected dofs as well. Defaults to False.
normalize_amplitude (Union[float, complex]) – If normalization is active, we can scale the overall magnitude of the eigenvector by this value. Defaults to 1.
normalize_max (bool) – If True, we normalize by the maximum magnitude of the listed dofs, otherwise by the average magnitude. Defaults to True.

Returns:

The processed eigenvectors, return it as result of the method process_eigenvectors().

run(endtime, timestep=None, *, outstep=None, numouts=None, out_initially=None, temporal_error=None, outstep_relative_to_zero=True, spatial_adapt=0, startstep=None, maxstep=None, newton_solver_tolerance=None, do_not_set_IC=False, globally_convergent_newton=False, max_newton_iterations=None, starttime=None, suppress_resolve_after_adapt=False, max_newton_to_increase_time_step=None)[source]

Run the problem for a specified duration, potential with output calls and temporal and/or spatial adaptivity. All time quantities must be given in dimensional units, e.g. second, if you use e.g. set_scaling() with e.g. temporal=1*second for a dimensional problem.

Parameters:

endtime (Union[Expression, int, float]) – The end time of the simulation.
timestep (Union[Expression, int, float, None]) – The time step size. If not specified, it will be determined automatically, e.g. by the outstep.
outstep (Union[Expression, int, float, None, bool]) – The time interval between outputs. If set to True, outputs will be generated at each time step. If set to False, no outputs will be generated. If not specified, it will be set to the value of timestep.
numouts (Optional[int]) – The number of outputs to generate. If specified, it will override the value of outstep.
out_initially (Optional[bool]) – Whether to generate an output at the initial time. If not specified, it will be set to True if outstep is not False, otherwise it will be set to False.
temporal_error (Optional[float]) – The temporal error tolerance. If specified, it will be used to control the time step size.
outstep_relative_to_zero (bool) – Whether the outstep is relative to the initial time or the current time. If set to True, the outstep will be relative to the initial time. If set to False, the outstep will be relative to the current time.
spatial_adapt (int) – The level of spatial adaptation. If specified, it will be used to control the spatial refinement level.
startstep (Union[Expression, int, float, None]) – The time step size at the start of the simulation (for temporal adaptivity). If specified, it will override the value of timestep.
maxstep (Union[Expression, int, float, None]) – The maximum time step size. If specified, it will be used to limit the time step size during temporal adaptivity.
newton_solver_tolerance (Optional[float]) – The tolerance for the Newton solver. If specified, it will be used to control the convergence of the solver during this run call.
do_not_set_IC (bool) – Whether to set the initial condition. If set to True, the initial condition will not be set.
globally_convergent_newton (bool) – Whether to use a globally convergent Newton solver. If set to True, a globally convergent Newton solver will be used.
max_newton_iterations (Optional[int]) – The maximum number of iterations for the Newton solver. If specified, it will override to limit the number of iterations for this run call.
starttime (Union[Expression, int, float, None]) – The start time of the simulation. If specified, it will override the current time.
suppress_resolve_after_adapt – Whether to suppress the resolve after adaptation. If set to True, the resolve after adaptation will be suppressed.
max_newton_to_increase_time_step (Optional[int]) – The maximum number of Newton iterations to increase the time step size. If specified, the adaptive time step will only be increased if the number of Newton iterations is less than this value.

Return type:

Union[Expression, int, float]

Returns:

The final time of the simulation after the run call.

Raises:

ValueError – If endtime is not specified.
ValueError – If outstep and numouts are specified simultaneously.
RuntimeError – If a suitable time step cannot be determined.

set_c_compiler(compiler_or_name)[source]

Selects the C compiler for the problem. “tcc” is fast in compilation, but slower in execution. Good for setting up a problem class. “system” is slower in compilation, but faster in execution. Good for running the final problem. set_c_compiler(“system”).optimize_for_max_speed() makes it even faster by using compiler flags for maximum speed.

Parameters:: compiler_or_name (Union[str, BaseCCompiler]) – The C compiler to use (“tcc” or “system” at the moment).
Returns:: The C compiler that was set.
Return type:: Union[_pyoomph.CCompiler, BaseCCompiler]

set_coordinate_system(csys)[source]

Set the default coordinate system at problem level. You can specify coordinate systems also at equation level, but if you don’t do, the coordinate system will default to this one.

Parameters:: csys (Union[Literal["axisymmetric","axisymmetric_flipped","cartesian","radialsymmetric"],BaseCoordinateSystem]) – The coordinate system to set as default.
Raises:: RuntimeError – Raised in case we do not set a valid coordinate system

set_current_time(val, dimensional=True, as_float=False)[source]

Set the current time of the problem.

Parameters:

val (ExpressionOrNum) – The value of the time to set.
dimensional (bool, optional) – Flag indicating whether the time is dimensional or not. Defaults to True.
as_float (bool, optional) – Flag indicating whether to convert the time to a float. Defaults to False.

Raises:

ValueError – If the nondimensional time is not a number.
ValueError – If the dimensional time cannot be nondimensionalized.

Returns:

None

set_eigensolver(solv)[source]

Set the eigensolver backend. “scipy”, “pardiso”, “slepc” are available (the latter two only if the packages MKL and/or petsc4py/slepc4py are installed)

Returns:: The eigenproblem solver instance after setting

set_initial_condition(ic_name='', all_unset_dofs_to_zero=False)[source]

Set the initial condition for the problem.

Parameters:

ic_name (str, optional) – Name of the initial condition. Multiple initial conditions can be defined in the problem definition by using the optional argument InitialConditions. Defaults to “”.
all_unset_dofs_to_zero (bool, optional) – Flag indicating whether to set all unset degrees of freedom, i.e. without any InitialCondition with this ic_name, to zero. Defaults to False.

set_linear_solver(solv)[source]

Set the linear solver backend. “scipy”, “umfpack”, “pardiso”, “petsc” are available (the latter two only if the packages MKL and/or petsc4py are installed)

Returns:: The linear solver instance after setting

set_output_directory(d)[source]

Change the output directory of the problem. Note: It should not be changed after the problem is initialised.

Parameters:: d (str) – Output directory
Return type:: None

set_scaling(**kwargs)[source]

Set the scaling factors for the problem variables for nondimensionalization. You can provide also scaling at equation level, but if not set there, it will ultimately default to the problem level scaling. Particular scales are "temporal" for the time and "spatial" for the spatial coordinates.

Parameters:

**kwargs (Union[Expression, int, float, str]) –

Keyword arguments specifying the scaling factors. The keys are the variable names, and the values are either numerical scaling factors or string expressions. In the latter case, we can set one scaling to another one, e.g.

set_scaling(u=1*meter/second,v="u")

would set the scaling of “v” to the one of “u”

Return type:

None

setup_for_stability_analysis(analytic_hessian=True, use_hessian_symmetry=True, improve_pitchfork_on_unstructured_mesh=False, improve_pitchfork_coordsys=None, improve_pitchfork_position_coordsys=None, shared_shapes_for_multi_assemble=None, azimuthal_stability=None, additional_cartesian_mode=None)[source]

Sets up the problem for stability analysis, e.g. for improved pitchfork tracking on unsymmetric meshes, azimuthal stability, etc. Arguments which are None are not changed.

Parameters:

analytic_hessian (bool, optional) – Flag indicating whether to use an analytically derived symbolical Hessian. Defaults to True.
use_hessian_symmetry (bool, optional) – Flag indicating whether to use symmetry in the Hessian. Defaults to True.
improve_pitchfork_on_unstructured_mesh (bool, optional) – Flag indicating whether to improve pitchfork tracking on unsymmetric meshes. Defaults to False.
improve_pitchfork_coordsys (OptionalCoordinateSystem, optional) – Coordinate system for improving pitchfork tracking unsymmetric meshes. Defaults to None.
improve_pitchfork_position_coordsys (OptionalCoordinateSystem, optional) – Coordinate system for improving pitchfork position space. Defaults to None.
shared_shapes_for_multi_assemble (Optional[bool], optional) – Flag indicating whether to use shared shapes for multi-assemble. Defaults to None.
azimuthal_stability (Optional[bool], optional) – Flag indicating whether to set up azimuthal stability code. Defaults to None.

solve(*, spatial_adapt=0, timestep=None, shift_values=True, temporal_error=None, max_newton_iterations=None, newton_relaxation_factor=None, suppress_resolve_after_adapt=False, newton_solver_tolerance=None, do_not_set_IC=False, globally_convergent_newton=False)[source]

Solves the problem stationary, unless a timestep is given. In that case, the time step is taken.

Parameters: - spatial_adapt (int): The level of spatial adaptation. Default is 0. - timestep (Union[ExpressionNumOrNone, List[ExpressionNumOrNone]]): The time step(s) for the transient solve. Can be a single value or a list of values. Default is None, meaning stationary solve without advancing in time. - shift_values (bool): Whether to shift the values during the solve, i.e. shifting the history value buffer. Default is True. - temporal_error (Optional[float]): The temporal error for adaptive time stepping. Default is None. - max_newton_iterations (Optional[int]): Override the maximum number of Newton iterations. Default is None. - newton_relaxation_factor (Optional[float]): Override the relaxation factor for the Newton solver. Default is None. - suppress_resolve_after_adapt (bool): Whether to suppress resolving after adaptation. Default is False. - newton_solver_tolerance (Optional[float]): Override the tolerance for the Newton solver. Default is None. - do_not_set_IC (bool): Whether to not set the initial condition in the first call. Default is False. - globally_convergent_newton (bool): Whether to use globally convergent Newton solver. Default is False.

Returns: - ExpressionOrNum: The current time after solving.

Return type:: Union[Expression, int, float]

solve_eigenproblem(n, shift=0, quiet=False, azimuthal_m=None, normal_mode_k=None, normal_mode_L=None, report_accuracy=False, sort=True, which='LM', OPpart=None, v0=None, filter=None, target=None)[source]

Solves the associated generalized eigenproblem for the given number of eigenvalues and eigenvectors.

Parameters:

n (int) – The number of eigenvalues and eigenvectors to compute.
shift (Union[float, complex, None], optional) – The shift applied for shift-inverted approaches to solve the eigenproblem. Defaults to 0.
quiet (bool, optional) – If True, suppresses the output. Defaults to False.
azimuthal_m (Optional[Union[int, List[int]]], optional) – The azimuthal mode number(s) for axial symmetry breaking. Defaults to None, i.e. the axisymmetric mode.
normal_mode_k (Union[Expression, int, float, List[Union[Expression, int, float]], None]) – The wave number(s) for an additional direction in Cartesian coordinates. Defaults to None, i.e. the base mode.
normal_mode_L (Union[Expression, int, float, List[Union[Expression, int, float]], None]) – The periodic length for an additional direction in Cartesian coordinates. Defaults to None, i.e. the base mode.
report_accuracy (bool, optional) – If True, reports the accuracy of the computed eigenvalues. Defaults to False.
sort (bool, optional) – If True, sorts the eigenvalues in ascending order. Defaults to True.
which ("EigenSolverWhich", optional) – The type of eigenvalues to compute. Defaults to “LM”.
OPpart (Optional[Literal["r", "i"]], optional) – The part of the operator to use. Defaults to None.
v0 (Optional[Union[NPFloatArray, NPComplexArray]], optional) – The initial guess for the eigenvectors. Defaults to None.
filter (Optional[Callable[[complex], bool]], optional) – A function to filter the computed eigenvalues. Only the eigenvalues for which the filter returns True will be kept. Defaults to None.
target (Optional[complex], optional) – The target eigenvalue. Defaults to None.

Returns:

A tuple containing the computed eigenvalues and eigenvectors.

Return type:

Tuple[NPComplexArray, NPComplexArray]

switch_to_hopf_orbit(eps=0.01, dparam=None, NT=30, mode='collocation', order=3, GL_order=-1, T_constraint='phase', amplitude_factor=1, FD_delta=1e-05, FD_param_delta=0.001, do_solve=True, solve_kwargs={}, check_collapse_to_stationary=True, orbit_amplitude=None, patch_number_of_nodes=True)[source]

After solving for a Hopf bifurcation by bifurcation tracking, this method will calculate the first Lyapunov exponent and initializes a good guess for the tracking of the periodic orbits originating at this Hopf bifurcation.

It is best to call it like:

with problem.switch_to_hopf_orbit(…) as orbit:
…

to deactivate orbit tracking after the with-statement.

Parameters:

eps (float) – A small number to construct the initial guess of the orbit and shift the parameter accordingly. Defaults to 0.01.
dparam (Optional[float]) – Optional parameter shift. If given and orbit_amplitude is also given, eps is ignored. Defaults to None.
NT (int) – Number of discrete time steps to consider for the orbit. Defaults to 30.
mode (Literal['collocation', 'floquet', 'central', 'BDF2', 'bspline']) – Selects the time discretization and interpolation mode. Defaults to “collocation”.
order (int) – Selects the order of the time discretization method. Defaults to 3.
GL_order (int) – Selects the Gauss-Legendre integration order for some time discretization modes. Defaults to -1, which is auto-select depending on the order.
T_constraint (Literal['phase', 'plane']) – Either use the “plane” or the “phase” constraint as equation for T. Defaults to “phase”.
amplitude_factor (float) – Additional multiplicative factor for the amplitude of the orbit guess. Defaults to 1.
FD_delta (float) – Finite difference step for the third order calculations used in the determination of the first Lyapunov coefficient. Defaults to 1e-5.
FD_param_delta – Finite difference step to determine the change of the real part of the eigenvalue with respect to the parameter. Defaults to 1e-3.
do_solve (bool) – Solve the orbit guess. Defaults to True.
solve_kwargs (Dict[str, Any]) – Additional keywords arguments to pass to the solve method for the initial solve. Defaults to {}.
check_collapse_to_stationary (bool) – Since an orbit can collapse to the stationary Hopf branch, we can check for it to make sure we are actually on an orbit. Defaults to True.
orbit_amplitude (Optional[float]) – Amplitude for the orbit. If set together with dparam, eps is ignored. Defaults to None.
patch_number_of_nodes (bool) – Depending on the order, we might have to slightly modify NT to have the right number of time nodes. Defaults to True.

Returns:

The periodic orbit object

Return type:

PeriodicOrbit

warn_about_unused_global_parameters: Union[bool, Literal['error']]: When set to True, we warn about unused global parameters for arclength continuation or bifurcation tracking. If set to “error”, we raise an error.

class pyoomph.generic.ResidualContribution(r, destination=None)[source]

Bases: BaseEquations

A class to add an arbitrary residual contribution to the equations. This is useful to add additional terms to the equations that are not covered by the standard weak formulation. Essentially, it just adds r to the residuals.

Parameters:

r (Union[Expression, int, float, str]) – The residual to add (can be e.g. a weak() contribution).
destination (Optional[str]) – The residual destination of the weak contribution. Can be used to define multiple residuals.

class pyoomph.generic.ScalarField(name, space, scale=None, testscale=None, residual=None)[source]

Bases: Equations

Introduces a scalar field with the given name and the given space. Residuals can be either added in the constructor or by combining with WeakContribution.

Parameters:

name (str) – Name of the scalar field
space (Literal['C1', 'C1TB', 'C2', 'C2TB', 'D1', 'D1TB', 'D2', 'D2TB', 'DL', 'D0']) – Space of the scalar field
scale (Union[Expression, int, float, None]) – Optional scaling of the field (default is 1)
testscale (Union[Expression, int, float, None]) – Optional scaling of the test function (default is 1)
residual (Union[Expression, int, float, None]) – Optional residual to be added. Formulate it in terms of the scalar field and the test function.

class pyoomph.generic.VectorField(name, space, scale=None, testscale=None, residual=None, dim=None)[source]

Bases: Equations

Introduces a vector field with the given name and the given space. Residuals can be either added in the constructor or by combining with WeakContribution.

Parameters:

name (str) – Name of the scalar field
space (Literal['C1', 'C1TB', 'C2', 'C2TB', 'D1', 'D1TB', 'D2', 'D2TB', 'DL', 'D0']) – Space of the scalar field
scale (Union[Expression, int, float, None]) – Optional scaling of the field (default is 1)
testscale (Union[Expression, int, float, None]) – Optional scaling of the test function (default is 1)
residual (Union[Expression, int, float, None]) – Optional residual to be added. Formulate it in terms of the scalar field and the test function.
dim (Optional[int]) – Vector dimension. If not set, it will be taken by the dimension of the mesh coordinates, i.e. the nodal dimension

class pyoomph.generic.WeakContribution(a, b, dimensional_dx=False, lagrangian=False, coordinate_system=None, destination=None)[source]

Bases: BaseEquations

A class to add an arbitrary weak contribution to the equations. This is useful to add additional terms to the equations that are not covered by the standard weak formulation. Essentially, it just adds weak(a,b) to the residuals.

Parameters:

a (Union[Expression, int, float, str]) – The lhs of the weak() contribution.
b (Union[Expression, str]) – The rhs (usually a testfunction()) of the weak() contribution.
dimensional_dx (bool) – If set to True, the weak contribution is treated as a dimensional contribution, i.e. spatial integration dx will carry dimension.
lagrangian (bool) – If set to True, the weak contribution is integrated in the Lagrangian frame of reference.
coordinate_system (Optional[BaseCoordinateSystem]) – The coordinate system in which the weak contribution is defined. If not set, the coordinate system of the equations or the problem is used.
destination (Optional[str]) – The residual destination of the weak contribution. Can be used to define multiple residuals.