Implementation of InverseIceVelocitySolver for the vertically integrated (SSA,SSA*,L1L2) stress models. More...

#include <InverseVerticallyIntegratedVelocitySolver.H>

Inheritance diagram for InverseVerticallyIntegratedVelocitySolver:

[legend]

Collaboration diagram for InverseVerticallyIntegratedVelocitySolver:

[legend]

Classes
class	Configuration
	Data read through ParmParse. More...

Public Member Functions
	InverseVerticallyIntegratedVelocitySolver ()

	~InverseVerticallyIntegratedVelocitySolver ()

void	setPreviousTime (Real a_time)

void	define (const ProblemDomain &a_coarseDomain, ConstitutiveRelation a_constRel, BasalFrictionRelation a_FrictionRel, const Vector< DisjointBoxLayout > &a_grids, const Vector< int > &a_refRatio, const RealVect &a_dxCrse, IceThicknessIBC *a_bc, int a_numLevels)
	(Re-)define IceVelocitySolver object for a given mesh and set of physics rules More...

void	define (const AmrIceBase &a_amrIce, const ProblemDomain &a_coarseDomain, ConstitutiveRelation a_constRel, BasalFrictionRelation a_FrictionRel, const Vector< DisjointBoxLayout > &a_grids, const Vector< int > &a_refRatio, const RealVect &a_dxCrse, IceThicknessIBC *a_bc, int a_numLevels)

int	solve (Vector< LevelData< FArrayBox > * > &a_horizontalVel, Vector< LevelData< FArrayBox > * > &a_calvedIce, Vector< LevelData< FArrayBox > * > &a_addedIce, Vector< LevelData< FArrayBox > * > &a_removedIce, Real &a_initialResidualNorm, Real &a_finalResidualNorm, const Real a_convergenceMetric, const Vector< LevelData< FArrayBox > * > &a_rhs, const Vector< LevelData< FArrayBox > * > &a_C, const Vector< LevelData< FArrayBox > * > &a_C0, const Vector< LevelData< FArrayBox > * > &a_A, const Vector< LevelData< FArrayBox > * > &a_muCoef, Vector< RefCountedPtr< LevelSigmaCS > > &a_coordSys, Real a_time, int a_lbase, int a_maxLevel)
	Compute the horizontal, cell centered ice velocity ( u(x,y), v(x,y) ) [VERTICALLY INTEGRATED MODELS!]. More...

void	computeObjectiveAndGradient (Real &a_fm, Real &a_fp, Vector< LevelData< FArrayBox > * > &a_g, const Vector< LevelData< FArrayBox > * > &a_x, bool a_inner)

void	restart ()
	carry out whatver operations are required prior to a restart More...

void	create (Vector< LevelData< FArrayBox > * > &a_a, const Vector< LevelData< FArrayBox > * > &a_b)

void	preCond (Vector< LevelData< FArrayBox > * > &a_s, const Vector< LevelData< FArrayBox > * > &a_r)

void	free (Vector< LevelData< FArrayBox > *> &a_data)
	need to make this explicit More...

void	plus (Vector< LevelData< FArrayBox > * > &a_a, Real a_r)

void	scale (Vector< LevelData< FArrayBox > * > &a_x, const Real a_s)

void	assign (Vector< LevelData< FArrayBox > * > &a_y, const Vector< LevelData< FArrayBox > * > &a_x)

void	incr (Vector< LevelData< FArrayBox > * > &y, const Vector< LevelData< FArrayBox > * > &x, Real s)

Real	dotProduct (Vector< LevelData< FArrayBox > * > &a_y, const Vector< LevelData< FArrayBox > * > &a_x)

int	nDoF (const Vector< LevelData< FArrayBox > * > &x)

BasalFriction *	basalFriction ()
	Create a BasalFriction object containing the solution . More...

MuCoefficient *	muCoefficient ()
	is the basal friction an improvement on the initial guess? More...

Public Member Functions inherited from IceVelocitySolver
	IceVelocitySolver ()

virtual	~IceVelocitySolver ()

virtual void	setVerbosity (int a_verbosity)

virtual void	setMaxIterations (int a_max_iter)
	if relevant, set max number of solver iterations More...

virtual void	setTolerance (Real a_tolerance)
	if relevant, set solver tolerance More...

virtual bool	isDefined ()

Public Member Functions inherited from ObjectiveWithGradient< Vector< LevelData< FArrayBox > *> >
virtual	~ObjectiveWithGradient ()
	virtual destructor More...

virtual Real	dotProduct (Vector< LevelData< FArrayBox > * > &a, const Vector< LevelData< FArrayBox > * > &b)=0

virtual void	computeObjectiveAndGradient (Real &fm, Real &fp, Vector< LevelData< FArrayBox > * > &g, const Vector< LevelData< FArrayBox > * > &x, bool inner)=0

virtual void	assign (Vector< LevelData< FArrayBox > * > &y, const Vector< LevelData< FArrayBox > * > &x)=0
	Copy x to y. More...

virtual void	incr (Vector< LevelData< FArrayBox > * > &y, const Vector< LevelData< FArrayBox > * > &x, Real s)=0
	Set y = y + s * x. More...

virtual void	create (Vector< LevelData< FArrayBox > * > &a, const Vector< LevelData< FArrayBox > * > &b)=0
	resize a / allocate storage to conform with b More...

virtual void	scale (Vector< LevelData< FArrayBox > * > &a, const Real s)=0
	set a = a * s More...

virtual void	preCond (Vector< LevelData< FArrayBox > * > &a, const Vector< LevelData< FArrayBox > * > &b)=0
	apply precoditioner K, so a = Kb More...

virtual int	nDoF (const Vector< LevelData< FArrayBox > * > &x)=0
	number of degrees-of-freedom (maximum number of iterations before restart) More...

Additional Inherited Members
Static Public Member Functions inherited from InverseIceVelocitySolver
static BasalFriction *	basalFriction (const Vector< LevelData< FArrayBox > *> &a_C, const Vector< int > &a_ratio, const RealVect &a_dxCrse)

static MuCoefficient *	muCoefficient (const Vector< LevelData< FArrayBox > *> &a_mu, const Vector< int > &a_ratio, const RealVect &a_dxCrse)

Protected Attributes inherited from IceVelocitySolver
bool	m_isDefined

int	m_verbosity

Detailed Description

Implementation of InverseIceVelocitySolver for the vertically integrated (SSA,SSA*,L1L2) stress models.

Computes a velocity $u$ , a basal friction coefficient $C$ and a stiffening coefficient $\phi$ which mimimises an objective function

$J = || f(u,u_o) ||^2_2 + || g(\nabla \cdot ( u h) , \nabla \cdot ( u_o h_o) ) ||^2_2 + \alpha _{0} || \nabla C ||^2_2 + \alpha _{1} || \nabla \phi ||^2_2 + \alpha _{2} || \nabla X_0 ||^2_2 + \alpha _{3} || \nabla X_1 ||^2_2 + \alpha _{4} || X_0 ||^2_2 + \alpha _{5} || X_1 ||^2_2$

This is essentially the inverse problem described in e.g Cornford et al, 2015, The Cryopshere.

Because $C$ and $\phi$ are positive, they are expressed as $C = C_0 \exp (X_0)$ and $\phi = \phi_0 \exp (X_1)$ .

The functions f,g are composed from a few parts to:

weight observations according to their accuracy. For example, f(u,u_o) can be zero in regions where u_o is unknown
choose a particular objective function, e.g, f (u,u_o) = |u| - |u_o| can be selected to match speeds rather than velocities.

The first part of J is minimal when the model and observed velocities match, given the form of f.

The second part is minimal when the model and observed flux divergence match. As this is often a noisy field, there is an option to smooth it

The remaining parts of $J$ are used to regularize the solution. Any or all the $\alpha _{n}$ can be zero. Typically, in time-independent problems only $\alpha _{0}$ and $\alpha _{1}$ are non-zero - this is Tikhonov regularization.

New to the code, and experimental, is some support for time independent problems. For now this is limited to recomputing $C(t)$ and $\phi(t)$ to minimise a $ J(t) $ . The factors $\alpha _{4} , \alpha _{5}$ can be made to depend on the time-step, so that $C(t + \Delta t)$ will tend to be closer to $C(t)$ as $\Delta t \rightarrow 0$ .

J is reduced in a series of nonlinear conjugate gradient (CG) iterations. To that end, the class implements the interface implied by the function template CGOptimize. The major method is computeObjectiveAndGradient, which calculates J and its directional derivatives with respect to X_0 and X_1. This function also writes the current state of the problem to hdf5 files, either every 'outer' CG iteration, or even more frequently, on every call, which will correspond to CG's 'inner' (secant) iterations

Uses JFNKSolver to solve a number of internal problems.

Constructor & Destructor Documentation

◆ InverseVerticallyIntegratedVelocitySolver()

InverseVerticallyIntegratedVelocitySolver::InverseVerticallyIntegratedVelocitySolver ( )

◆ ~InverseVerticallyIntegratedVelocitySolver()

InverseVerticallyIntegratedVelocitySolver::~InverseVerticallyIntegratedVelocitySolver ( )

Member Function Documentation

◆ assign()

void InverseVerticallyIntegratedVelocitySolver::assign	(	Vector< LevelData< FArrayBox > * > &	a_y,
		const Vector< LevelData< FArrayBox > * > &	a_x
	)

References incr().

Referenced by computeObjectiveAndGradient(), dotProduct(), preCond(), scale(), and solve().

◆ basalFriction()

BasalFriction* InverseVerticallyIntegratedVelocitySolver::basalFriction ( )

inlinevirtual

Create a BasalFriction object containing the solution $C$ .

Implements InverseIceVelocitySolver.

References InverseIceVelocitySolver::basalFriction().

◆ computeObjectiveAndGradient()

void InverseVerticallyIntegratedVelocitySolver::computeObjectiveAndGradient	(	Real &	a_fm,
		Real &	a_fp,
		Vector< LevelData< FArrayBox > * > &	a_g,
		const Vector< LevelData< FArrayBox > * > &	a_x,
		bool	a_inner
	)

Referenced by nDoF().

◆ create()

void InverseVerticallyIntegratedVelocitySolver::create	(	Vector< LevelData< FArrayBox > * > &	a_a,
		const Vector< LevelData< FArrayBox > * > &	a_b
	)

References preCond().

◆ define() [1/2]

void InverseVerticallyIntegratedVelocitySolver::define	(	const ProblemDomain &	a_coarseDomain,
		ConstitutiveRelation *	a_constRel,
		BasalFrictionRelation *	a_FrictionRel,
		const Vector< DisjointBoxLayout > &	a_vectGrids,
		const Vector< int > &	a_vectRefRatio,
		const RealVect &	a_dxCrse,
		IceThicknessIBC *	a_bc,
		int	a_numLevels
	)

inlinevirtual

(Re-)define IceVelocitySolver object for a given mesh and set of physics rules

Parameters

a_coarseDomain	domain boundaries on the coarsest level
a_constRel	physics rule (e.g Glen's flow law) relating englacial stress and strian-rate
a_FrictionRel	physics rule (e.g Weertman law) relating basal stress and strain-rate
a_vectGrids	mesh layout for each level
a_vectGrids	refinement ratio for each level
a_dxCrse	mesh spacing on coarsest level
a_bc	ice thickness and velocity boundary conditions
a_numLevels	number of mesh levels.

Implements IceVelocitySolver.

References IceVelocitySolver::define(), and IceVelocitySolver::solve().

Referenced by AmrIce::defineSolver(), and InverseVerticallyIntegratedVelocitySolver::Configuration::parse().

◆ define() [2/2]

void InverseVerticallyIntegratedVelocitySolver::define	(	const AmrIceBase &	a_amrIce,
		const ProblemDomain &	a_coarseDomain,
		ConstitutiveRelation *	a_constRel,
		BasalFrictionRelation *	a_FrictionRel,
		const Vector< DisjointBoxLayout > &	a_grids,
		const Vector< int > &	a_refRatio,
		const RealVect &	a_dxCrse,
		IceThicknessIBC *	a_bc,
		int	a_numLevels
	)

References InverseVerticallyIntegratedVelocitySolver::Configuration::parse(), and solve().

◆ dotProduct()

Real InverseVerticallyIntegratedVelocitySolver::dotProduct	(	Vector< LevelData< FArrayBox > * > &	a_y,
		const Vector< LevelData< FArrayBox > * > &	a_x
	)

References assign(), and scale().

Referenced by incr().

◆ free()

void InverseVerticallyIntegratedVelocitySolver::free ( Vector< LevelData< FArrayBox > *> & a_data )

inlinevirtual

need to make this explicit

Implements ObjectiveWithGradient< Vector< LevelData< FArrayBox > *> >.

◆ incr()

void InverseVerticallyIntegratedVelocitySolver::incr	(	Vector< LevelData< FArrayBox > * > &	y,
		const Vector< LevelData< FArrayBox > * > &	x,
		Real	s
	)

References dotProduct().

Referenced by assign().

◆ muCoefficient()

MuCoefficient* InverseVerticallyIntegratedVelocitySolver::muCoefficient ( )

inlinevirtual

is the basal friction an improvement on the initial guess?

Create a MuCoefficient object containing the solution $\phi$

Implements InverseIceVelocitySolver.

References InverseIceVelocitySolver::muCoefficient().

◆ nDoF()

int InverseVerticallyIntegratedVelocitySolver::nDoF ( const Vector< LevelData< FArrayBox > * > & x )

◆ plus()

void InverseVerticallyIntegratedVelocitySolver::plus	(	Vector< LevelData< FArrayBox > * > &	a_a,
		Real	a_r
	)

References scale().

Referenced by computeObjectiveAndGradient(), preCond(), and solve().

◆ preCond()

void InverseVerticallyIntegratedVelocitySolver::preCond	(	Vector< LevelData< FArrayBox > * > &	a_s,
		const Vector< LevelData< FArrayBox > * > &	a_r
	)

References assign(), and plus().

Referenced by create().

◆ restart()

void InverseVerticallyIntegratedVelocitySolver::restart ( )

inlinevirtual

carry out whatver operations are required prior to a restart

Implements ObjectiveWithGradient< Vector< LevelData< FArrayBox > *> >.

◆ scale()

void InverseVerticallyIntegratedVelocitySolver::scale	(	Vector< LevelData< FArrayBox > * > &	a_x,
		const Real	a_s
	)

References assign().

Referenced by computeObjectiveAndGradient(), dotProduct(), and plus().

◆ setPreviousTime()

void InverseVerticallyIntegratedVelocitySolver::setPreviousTime ( Real a_time )

inlinevirtual

Implements InverseIceVelocitySolver.

◆ solve()

int InverseVerticallyIntegratedVelocitySolver::solve	(	Vector< LevelData< FArrayBox > * > &	a_horizontalVel,
		Vector< LevelData< FArrayBox > * > &	a_calvedIce,
		Vector< LevelData< FArrayBox > * > &	a_addedIce,
		Vector< LevelData< FArrayBox > * > &	a_removedIce,
		Real &	a_initialResidualNorm,
		Real &	a_finalResidualNorm,
		const Real	a_convergenceMetric,
		const Vector< LevelData< FArrayBox > * > &	a_rhs,
		const Vector< LevelData< FArrayBox > * > &	a_C,
		const Vector< LevelData< FArrayBox > * > &	a_C0,
		const Vector< LevelData< FArrayBox > * > &	a_A,
		const Vector< LevelData< FArrayBox > * > &	a_muCoef,
		Vector< RefCountedPtr< LevelSigmaCS > > &	a_coordSys,
		Real	a_time,
		int	a_lbase,
		int	a_maxLevel
	)

virtual

Compute the horizontal, cell centered ice velocity ( u(x,y), v(x,y) ) [VERTICALLY INTEGRATED MODELS!].

Parameters

a_horizontalVel	on input an initial guess ( u_0(x,y), v_0(x,y) ), on output ( u(x,y), v(x,y) )
a_initialResidualNorm	residual norm given the initial guess
a_finalResidualNorm	residual norm at the end of the calculation
a_convergenceMetric	desired a_finalResidualNorm/a_initialResidualNorm
a_rhs	right hand side of the stress balance equation: the driving stress (x,y) components
a_C	main basal friction coefficient : $\tau _b = f_b(C,u) + C_0 u$
a_C0	auxilliary basal friction coefficient : $\tau _b =f_b(C,u) + C_0 u$
a_A	flow rate factor: one component per layer
a_muCoef	stiffness factor: one component
a_coordSys	ice sheet geoemtry
a_time	current time
a_lbase	lowest mesh level (probably doesn't work!)
a_maxLevel	highest mesh level

Returns: 0 if solver converges successfully

Implements IceVelocitySolver.

Referenced by define().

The documentation for this class was generated from the following files:

Classes

Public Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ InverseVerticallyIntegratedVelocitySolver()

◆ ~InverseVerticallyIntegratedVelocitySolver()

Member Function Documentation

◆ assign()

◆ basalFriction()

◆ computeObjectiveAndGradient()

◆ create()

◆ define() [1/2]

◆ define() [2/2]

◆ dotProduct()

◆ free()

◆ incr()

◆ muCoefficient()

◆ nDoF()

◆ plus()

◆ preCond()

◆ restart()

◆ scale()

◆ setPreviousTime()

◆ solve()