Inversion

Contents

rhea_inversion.h
rhea_inversion_obs_stress.h
rhea_inversion_obs_velocity.h
rhea_inversion_obs_viscosity.h
rhea_inversion_param.h
rhea_inversion_qoi.h

---

rhea_inversion.h [source]

Typedefs

typedef struct rhea_inversion_problem rhea_inversion_problem_t

Enums

enum rhea_inversion_krylov_solver_idx_t

RHEA_INVERSION

Inversion for Stokes parameters.

Values:

enumerator RHEA_INVERSION_KRYLOV_SOLVER_IDX_FORWARD = 0

enumerator RHEA_INVERSION_KRYLOV_SOLVER_IDX_ADJOINT

enumerator RHEA_INVERSION_KRYLOV_SOLVER_IDX_INCR_FORWARD

enumerator RHEA_INVERSION_KRYLOV_SOLVER_IDX_INCR_ADJOINT

enumerator RHEA_INVERSION_KRYLOV_SOLVER_N

Functions

void rhea_inversion_add_options(ymir_options_t *opt_sup)

Defines options and adds them as sub-options.

void rhea_inversion_perfmon_init(const int activate, const int skip_if_active)

Initializes performance counters.

void rhea_inversion_perfmon_print(sc_MPI_Comm mpicomm, const int print_wtime, const int print_n_calls, const int print_flops)

Prints statistics collected by performance monitors.

rhea_inversion_problem_t *rhea_inversion_new(rhea_stokes_problem_t *stokes_problem)

Creates/destroys an inverse problem.

void rhea_inversion_destroy(rhea_inversion_problem_t *inv_problem)

void rhea_inversion_set_txt_output(rhea_inversion_problem_t *inv_problem, char *txt_path)

Sets paths for writing output files.

void rhea_inversion_set_vtk_output(rhea_inversion_problem_t *inv_problem, char *vtk_path_vol, char *vtk_path_surf)

void rhea_inversion_solve(rhea_inversion_problem_t *inv_problem, const int use_initial_guess, ymir_vec_t *inv_parameter_vec)

Solves an inverse problem.

void rhea_inversion_solve_with_vel_obs(rhea_inversion_problem_t *inv_problem, const int use_initial_guess, ymir_vec_t *inv_parameter_vec, ymir_vec_t *vel_obs_surf, ymir_vec_t *vel_obs_weight_surf, const double vel_obs_add_noise_stddev)

Solves an inverse problem with given velocity observations at the surface.

---

rhea_inversion_obs_stress.h [source]

Enums

enum rhea_inversion_obs_stress_t

Values:

enumerator RHEA_INVERSION_OBS_STRESS_NONE = -1

enumerator RHEA_INVERSION_OBS_STRESS_VOLUME = 0

enumerator RHEA_INVERSION_OBS_STRESS_QOI_PLATE_BOUNDARY_NORMAL = 100

enumerator RHEA_INVERSION_OBS_STRESS_QOI_PLATE_BOUNDARY_TANGENTIAL_X = 101

enumerator RHEA_INVERSION_OBS_STRESS_QOI_PLATE_BOUNDARY_TANGENTIAL_Y = 102

enumerator RHEA_INVERSION_OBS_STRESS_QOI_PLATE_BOUNDARY_TANGENTIAL_Z = 103

Functions

void rhea_inversion_obs_stress_diff(ymir_vec_t *misfit_stress, ymir_vec_t *forward_vel, ymir_vec_t *obs_stress, ymir_vec_t *weight, rhea_stokes_problem_t *stokes_problem)

Compute the difference of the stresses: weight * (ObsOp(stress(vel,press)) - stress_obs) where ObsOp(stress(vel,press)) = 2 * visc * strainrate(vel) - press * I

double rhea_inversion_obs_stress_misfit(ymir_vec_t *forward_vel, ymir_vec_t *obs_stress, ymir_vec_t *weight, const rhea_inversion_obs_stress_t obs_type, rhea_stokes_problem_t *stokes_problem)

Computes misfit term (squared norm of the difference) of the of the difference in stresses.

double rhea_inversion_obs_stress_misfit_param_derivative(ymir_vec_t *forward_vel_press, ymir_vec_t *obs_stress, ymir_vec_t *weight, const rhea_inversion_obs_stress_t obs_type, rhea_stokes_problem_t *stokes_problem, ymir_stress_op_t *stress_op_param_derivative)

Computes derivative of the misfit term with respect to parameters.

void rhea_inversion_obs_stress_add_adjoint_rhs(ymir_vec_t *rhs_vel_mass, ymir_vec_t *forward_vel, ymir_vec_t *obs_stress, ymir_vec_t *weight, const rhea_inversion_obs_stress_t obs_type, rhea_stokes_problem_t *stokes_problem)

Adds contribution to the right-hand side for adjoint equations.

---

rhea_inversion_obs_velocity.h [source]

Enums

enum rhea_inversion_obs_velocity_t

Values:

enumerator RHEA_INVERSION_OBS_VELOCITY_NONE = -1

enumerator RHEA_INVERSION_OBS_VELOCITY_NORMAL = 0

enumerator RHEA_INVERSION_OBS_VELOCITY_TANGENTIAL

enumerator RHEA_INVERSION_OBS_VELOCITY_TANGENTIAL_ROTFREE

enumerator RHEA_INVERSION_OBS_VELOCITY_ALL

enumerator RHEA_INVERSION_OBS_VELOCITY_ALL_ROTFREE

enum rhea_inversion_obs_velocity_weight_t

Values:

enumerator RHEA_INVERSION_OBS_VELOCITY_WEIGHT_VALUES = 0

enumerator RHEA_INVERSION_OBS_VELOCITY_WEIGHT_INV_AREA_SQRT

enumerator RHEA_INVERSION_OBS_VELOCITY_WEIGHT_INV_AREA_LOG

enumerator RHEA_INVERSION_OBS_VELOCITY_WEIGHT_INV_AREA_LIN

Functions

void rhea_inversion_obs_velocity_generate(ymir_vec_t *vel_obs_surf, ymir_vec_t *weight_surf, const rhea_inversion_obs_velocity_t obs_type, const rhea_inversion_obs_velocity_weight_t weight_type, const double *weight_values, rhea_plate_options_t *plate_options, double *calculated_weight_values)

Generates the observational data.

void rhea_inversion_obs_velocity_diff(ymir_vec_t *misfit_surf, ymir_vec_t *vel_fwd_vol, ymir_vec_t *vel_obs_surf, ymir_vec_t *weight_surf, const rhea_inversion_obs_velocity_t obs_type, rhea_domain_options_t *domain_options)

Compute the difference of the velocities at the surface: weight * (ObsOp(vel) - vel_obs)

Compute the difference of the velocities at the surface: weight * (ObsOp(vel) - vel_obs) where ObsOp(vel) = ObsOp(vel_vol) = ProjectToSurface(vel_vol) = vel_surf

double rhea_inversion_obs_velocity_misfit(ymir_vec_t *vel_fwd_vol, ymir_vec_t *vel_obs_surf, ymir_vec_t *weight_surf, const rhea_inversion_obs_velocity_t obs_type, rhea_domain_options_t *domain_options)

Computes misfit term (squared norm of the difference) of the of velocities at the surface.

double rhea_inversion_obs_velocity_misfit_param_derivative(const rhea_inversion_obs_velocity_t obs_type)

Computes derivative of the misfit term with respect to parameters.

void rhea_inversion_obs_velocity_add_adjoint_rhs(ymir_vec_t *rhs_vel_mass, ymir_vec_t *vel_fwd_vol, ymir_vec_t *vel_obs_surf, ymir_vec_t *weight_surf, const rhea_inversion_obs_velocity_t obs_type, rhea_domain_options_t *domain_options)

Adds contribution to the right-hand side for adjoint equations.

void rhea_inversion_obs_velocity_incremental_adjoint_rhs(ymir_vec_t *rhs_vel_mass, ymir_vec_t *vel_incrfwd_vol, ymir_vec_t *weight_surf, const rhea_inversion_obs_velocity_t obs_type, rhea_domain_options_t *domain_options)

Computes the right-hand side (component) for incremental adjoint equations.

void rhea_inversion_obs_velocity_remove_normal(ymir_vec_t *vel_surf)

Removes normal/tangential components of a surface velocity field.

void rhea_inversion_obs_velocity_remove_tangential(ymir_vec_t *vel_surf)

---

rhea_inversion_obs_viscosity.h [source]

Enums

enum rhea_inversion_obs_viscosity_t

Values:

enumerator RHEA_INVERSION_OBS_VISCOSITY_NONE = -1

enumerator RHEA_INVERSION_OBS_VISCOSITY_AVERAGE_REGION = 0

enumerator RHEA_INVERSION_OBS_VISCOSITY_AVERAGE_UNDER_PLATES

Functions

rhea_domain_subset_column_t **rhea_inversion_obs_viscosity_new(int *n_columns, double **value, double **weight, const rhea_inversion_obs_viscosity_t obs_type, char *value_Pas_list, char *stddev_rel_list, rhea_plate_options_t *plate_options, rhea_viscosity_options_t *visc_options)

Creates observational data of avarage viscosities.

void rhea_inversion_obs_viscosity_destroy(rhea_domain_subset_column_t **column, const int n_columns, double *value, double *weight)

Destroys observational data of avarage viscosities.

double rhea_inversion_obs_viscosity_misfit(ymir_vec_t *forward_vel, const rhea_inversion_obs_viscosity_t obs_type, const int n_columns, rhea_domain_subset_column_t **column, const double *obs_val, const double *weight, rhea_stokes_problem_t *stokes_problem)

Computes the misfit term of average viscosities.

double rhea_inversion_obs_viscosity_misfit_param_derivative(const rhea_inversion_obs_viscosity_t obs_type)

Computes derivative of the misfit term with respect to parameters.

void rhea_inversion_obs_viscosity_add_adjoint_rhs(ymir_vec_t *rhs_vel_mass, ymir_vec_t *forward_vel, const rhea_inversion_obs_viscosity_t obs_type, const int n_columns, rhea_domain_subset_column_t **column, const double *obs_val, const double *weight, rhea_stokes_problem_t *stokes_problem)

Adds contribution to the right-hand side for adjoint equations.

---

rhea_inversion_param.h [source]

Typedefs

typedef struct rhea_inversion_param_options rhea_inversion_param_options_t

RHEA_INVERSION_PARAM

Parameters of inverse problems for mantle convection.

typedef struct rhea_inversion_param rhea_inversion_param_t

typedef double (*rhea_inversion_param_derivative_fn_t)(ymir_vec_t *forward_vel_press, ymir_vec_t *adjoint_vel_press, const int parameter_idx, const int derivative_type, const rhea_weakzone_label_t weak_label, ymir_vec_t *coeff_param_derivative, ymir_stress_op_t *stress_op_param_derivative, void *data)

Callback function that computes the derivative with respect to a single parameter.

Param parameter_idx:

[in] Index of parameter.

Enums

enum rhea_inversion_param_verbosity_t

Values:

enumerator RHEA_INVERSION_PARAM_VERBOSE_NONE

enumerator RHEA_INVERSION_PARAM_VERBOSE_REAL

enumerator RHEA_INVERSION_PARAM_VERBOSE_REAL_NONDIM

enumerator RHEA_INVERSION_PARAM_VERBOSE_REAL_NONDIM_DIM

Functions

void rhea_inversion_param_add_options(rhea_inversion_param_options_t *inv_param_options, ymir_options_t *opt_sup)

Defines options and adds them as sub-options.

rhea_inversion_param_t *rhea_inversion_param_new(rhea_stokes_problem_t *stokes_problem, rhea_inversion_param_options_t *inv_param_options)

Creates a new set of inversion parameters.

void rhea_inversion_param_destroy(rhea_inversion_param_t *inv_param)

Destroys a set of inversion parameters.

void rhea_inversion_param_print(ymir_vec_t *parameter_vec, const rhea_inversion_param_verbosity_t verbosity, rhea_inversion_param_t *inv_param)

Prints the active parameters with specified verbosity.

void rhea_inversion_param_pull_from_model(ymir_vec_t *parameter_vec, rhea_inversion_param_t *inv_param)

Sets the inversion parameters to the values from a model.

void rhea_inversion_param_push_to_model(ymir_vec_t *parameter_vec, rhea_inversion_param_t *inv_param)

Overwrites model parameters with values from the parameter inversion.

void rhea_inversion_param_set_initial_guess(ymir_vec_t *parameter_vec, rhea_inversion_param_t *inv_param, rhea_inversion_param_options_t *opt)

Sets the inversion parameters to initial values.

double rhea_inversion_param_derivative_min_max(const double model_val)

Computes the derivative of the parametrization function.

(d/dp *_convert_to_model_*) (*_convert_from_model_* (m))

where *_min_max () Used for lower and upper bounds of viscosity. *_scal () Used for scaling factors (pre-factors) of viscosity. *_arrh () Used for Arrhenius activation energy. *_n () Used for stress exponent. *_yield () Used for yield strength. *_weak () Used for weak zone factors.

double rhea_inversion_param_derivative_scal(const double model_val)

double rhea_inversion_param_derivative_arrh(const double model_val)

double rhea_inversion_param_derivative_n(const double model_val)

double rhea_inversion_param_derivative_yield(const double model_val)

double rhea_inversion_param_derivative_weak(const double model_val)

double rhea_inversion_param_prior(ymir_vec_t *parameter_vec, rhea_inversion_param_t *inv_param)

Computes the (squared norm of the) prior term for the objective functional.

void rhea_inversion_param_compute_gradient(ymir_vec_t *gradient_vec, ymir_vec_t *parameter_vec, ymir_vec_t *forward_vel_press, ymir_vec_t *adjoint_vel_press, const double prior_weight, rhea_inversion_param_t *inv_param, ymir_vec_t *gradient_adjoint_comp_vec, ymir_vec_t *gradient_prior_comp_vec, ymir_vec_t *gradient_callback_comp_vec, rhea_inversion_param_derivative_fn_t derivative_fn, void *derivative_fn_data)

Computes the gradient vector of the Stokes model w.r.t. model parameters.

double rhea_inversion_param_compute_gradient_norm(ymir_vec_t *gradient_vec, rhea_inversion_param_t *inv_param)

Computes the norm of the provided gradient vector.

void rhea_inversion_param_incremental_forward_rhs(ymir_vec_t *rhs_vel_mass, ymir_vec_t *gradient_direction, ymir_vec_t *forward_vel_press, rhea_inversion_param_t *inv_param)

Computes the right-hand side for incremental forward equations.

void rhea_inversion_param_apply_hessian(ymir_vec_t *param_vec_out, ymir_vec_t *param_vec_in, ymir_vec_t *forward_vel_press, ymir_vec_t *adjoint_vel_press, ymir_vec_t *incr_forward_vel_press, ymir_vec_t *incr_adjoint_vel_press, const double prior_weight, rhea_inversion_param_t *inv_param)

Applies the Hessian of the inverse problem to a parameter vector.

void rhea_inversion_param_prior_inv_cov(ymir_vec_t *inverse_covariance, rhea_inversion_param_t *inv_param)

Computes the (diagonal) inverse covariance of the prior.

int rhea_inversion_param_restrict_to_feasible(ymir_vec_t *parameter_vec, const double restrict_to_prior_stddev, rhea_inversion_param_t *inv_param)

Restricts the parameters to their feasible values (elementwise).

int rhea_inversion_param_restrict_step_length_to_feasible(ymir_vec_t *step_vec, ymir_vec_t *parameter_vec, const double restrict_to_prior_stddev, rhea_inversion_param_t *inv_param, double *step_length_new)

Modifies the step such that the new parameter vector, parameter_vec + step_vec, lies within the feasible set associated with the inversion parameters.

ymir_vec_t *rhea_inversion_param_vec_new(rhea_inversion_param_t *inv_param)

Creates/destroys a parameter vector.

ymir_vec_t *rhea_inversion_param_vec_new_perturb(rhea_inversion_param_t *inv_param, const double perturb_stddev)

void rhea_inversion_param_vec_destroy(ymir_vec_t *vec)

int rhea_inversion_param_vec_check_type(ymir_vec_t *vec, rhea_inversion_param_t *inv_param)

Checks whether a vector is of the right type.

int rhea_inversion_param_vec_is_valid(ymir_vec_t *vec, rhea_inversion_param_t *inv_param)

Checks entries of a vector.

void rhea_inversion_param_vec_print(ymir_vec_t *vec, rhea_inversion_param_t *inv_param)

Prints the active values of a parameter vector.

sc_dmatrix_t *rhea_inversion_param_vec_reduced_new(ymir_vec_t *vec, rhea_inversion_param_t *inv_param)

Creates/destroys a reduced parameter vector, where only the active entries are copied from the full parameter vector in the input argument.

void rhea_inversion_param_vec_reduced_destroy(sc_dmatrix_t *vec_reduced)

void rhea_inversion_param_vec_reduced_copy(ymir_vec_t *vec, sc_dmatrix_t *vec_reduced, rhea_inversion_param_t *inv_param)

Copy entries from a reduced to a full parameter vector.

void rhea_inversion_param_vec_reduced_apply_matrix(ymir_vec_t *out, sc_dmatrix_t *matrix_reduced, ymir_vec_t *in, rhea_inversion_param_t *inv_param)

Applies a reduced matrix to a full parameter vector.

int rhea_inversion_param_vec_reduced_solve_matrix(ymir_vec_t *sol, sc_dmatrix_t *matrix_reduced, ymir_vec_t *rhs, rhea_inversion_param_t *inv_param, int *n_iterations)

Solves a linear system with a reduced matrix and a right-hand side given by a full parameter vector.

double rhea_inversion_param_vec_reduced_inner_product(ymir_vec_t *vecL, ymir_vec_t *vecR, ymir_vec_t *weight, rhea_inversion_param_t *inv_param, const int normalize_wrt_size)

Computes inner product of two vectors.

double rhea_inversion_param_vec_reduced_ip(sc_dmatrix_t *vecL_reduced, sc_dmatrix_t *vecR_reduced, sc_dmatrix_t *weight_reduced, const int normalize_wrt_size)

int rhea_inversion_param_get_n_parameters(rhea_inversion_param_t *inv_param)

Returns the number of parameters.

int rhea_inversion_param_get_n_active(rhea_inversion_param_t *inv_param)

Returns the number of active parameters.

int *rhea_inversion_param_get_active(rhea_inversion_param_t *inv_param)

Gets the pointer to the activation mask.

struct rhea_inversion_param_options

#include <rhea_inversion_param.h>

RHEA_INVERSION_PARAM

Parameters of inverse problems for mantle convection.

Public Members

int min_a

int max_a

int upper_mantle_scaling_a

int upper_mantle_arrhenius_activation_energy_a

int lower_mantle_scaling_a

int lower_mantle_arrhenius_activation_energy_a

int stress_exponent_a

int yield_strength_a

int thickness_a

int thickness_class_slab_a

int thickness_class_ridge_a

int thickness_class_fracture_a

int thickness_const_a

int thickness_const_class_slab_a

int thickness_const_class_ridge_a

int thickness_const_class_fracture_a

int weak_factor_interior_a

int weak_factor_interior_class_slab_a

int weak_factor_interior_class_ridge_a

int weak_factor_interior_class_fracture_a

int weak_factor_interior_label_slab_a

int weak_factor_interior_label_ridge_a

int weak_factor_interior_label_fracture_a

char *weak_factor_interior_label_file_path_txt_a

double prior_mean_perturb_stddev

double prior_stddev_min

double prior_stddev_max

double prior_stddev_upper_mantle_scaling

double prior_stddev_upper_mantle_arrhenius_activation_energy

double prior_stddev_lower_mantle_scaling

double prior_stddev_lower_mantle_arrhenius_activation_energy

double prior_stddev_stress_exponent

double prior_stddev_yield_strength

double prior_stddev_thickness

double prior_stddev_thickness_const

double prior_stddev_weak_factor_interior

char *initial_guess_file_path_txt

double initial_guess_perturb_stddev

double initial_guess_shift_by_prior_stddev

---

rhea_inversion_qoi.h [source]

Typedefs

typedef struct rhea_inversion_qoi_options rhea_inversion_qoi_options_t

typedef struct rhea_inversion_qoi rhea_inversion_qoi_t

Enums

enum rhea_inversion_qoi_class_t

RHEA_INVERSION_QOI

Quantities of interest (QOI) of inverse problems for mantle convection.

Values:

enumerator RHEA_INVERSION_QOI_STRESS

Functions

void rhea_inversion_qoi_add_options(rhea_inversion_qoi_options_t *inv_qoi_options, ymir_options_t *opt_sup)

Defines options and adds them as sub-options.

rhea_inversion_qoi_t *rhea_inversion_qoi_new(rhea_stokes_problem_t *stokes_problem, rhea_inversion_qoi_options_t *inv_qoi_options)

Creates a new set of quantities of interest.

void rhea_inversion_qoi_destroy(rhea_inversion_qoi_t *inv_qoi)

Destroys a set of quantities of interest.

int rhea_inversion_qoi_exists(rhea_inversion_qoi_t *inv_qoi)

Checks if any QOI are set.

void rhea_inversion_qoi_stress_create_obs(rhea_inversion_obs_stress_t *stress_obs_type, ymir_vec_t **stress_obs_weight, rhea_inversion_qoi_t *inv_qoi, const int reduced_index)

void rhea_inversion_qoi_stress_clear_obs(ymir_vec_t *stress_obs_weight, rhea_inversion_qoi_t *inv_qoi)

sc_dmatrix_t *rhea_inversion_qoi_stress_vec_reduced_new(rhea_inversion_qoi_t *inv_qoi)

Creates/destroys a reduced QOI vector.

void rhea_inversion_qoi_stress_vec_reduced_destroy(sc_dmatrix_t *vec_reduced)

int rhea_inversion_qoi_stress_vec_set_entry(sc_dmatrix_t *vec_reduced, const double val, rhea_inversion_qoi_t *inv_qoi, const int stress_qoi_index, const int weakzone_label)

Sets an entry of a QOI vector. Returns index of which entry was set.

sc_dmatrix_t *rhea_inversion_qoi_stress_to_param_matrix_new(const int n_parameters, rhea_inversion_qoi_t *inv_qoi)

Creates/destroys a matrix mapping a reduced QOI vector to a reduced parameter space.

void rhea_inversion_qoi_stress_to_param_matrix_destroy(sc_dmatrix_t *mat)

int rhea_inversion_qoi_stress_to_param_matrix_set_column(sc_dmatrix_t *mat, sc_dmatrix_t *param_vec_reduced, rhea_inversion_qoi_t *inv_qoi, const int stress_qoi_index, const int weakzone_label)

Sets a column of a QOI-to-parameter matrix. Returns index of which column was set.

int rhea_inversion_qoi_get_count_all(rhea_inversion_qoi_t *inv_qoi)

Returns the number of QOI.

struct rhea_inversion_qoi_options

Public Members

char *stress_type_list

char *stress_activate_weakzone_label_file_path_txt