src/base/stressbalance/SSA.hh

Go to the documentation of this file.

// Copyright (C) 2004--2011 Jed Brown, Ed Bueler and Constantine Khroulev
//
// This file is part of PISM.
//
// PISM is free software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the Free Software
// Foundation; either version 2 of the License, or (at your option) any later
// version.
//
// PISM is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
// details.
//
// You should have received a copy of the GNU General Public License
// along with PISM; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

#ifndef _SSA_H_
#define _SSA_H_

#include "ShallowStressBalance.hh"
#include "PISMDiagnostic.hh"


class SSAStrengthExtension {
public:
  SSAStrengthExtension(const NCConfigVariable &c) : config(c) {
    min_thickness = config.get("min_thickness_strength_extension_ssa");
    constant_nu = config.get("constant_nu_strength_extension_ssa");
  }

  virtual ~SSAStrengthExtension() {}


  virtual PetscErrorCode set_notional_strength(PetscReal my_nuH) {
     if (my_nuH <= 0.0) SETERRQ(1,"nuH must be positive");
     constant_nu = my_nuH / min_thickness;
     return 0;
  }


  virtual PetscErrorCode set_min_thickness(PetscReal my_min_thickness) {
     if (my_min_thickness <= 0.0) SETERRQ(1,"min_thickness must be positive");
     PetscReal nuH = constant_nu * min_thickness;
     min_thickness = my_min_thickness;
     constant_nu = nuH / min_thickness;
     return 0;
  }

  virtual PetscReal get_notional_strength() const {
    return constant_nu * min_thickness;
  }

  virtual PetscReal get_min_thickness() const { return min_thickness; }

private:
  const NCConfigVariable &config;
  PetscReal  min_thickness, constant_nu;
};


class SSA;
typedef SSA * (*SSAFactory)(IceGrid &, IceBasalResistancePlasticLaw &, 
              IceFlowLaw &, EnthalpyConverter &, const NCConfigVariable &);



class SSA : public ShallowStressBalance
{
  friend class SSA_taud;
public:
  SSA(IceGrid &g, IceBasalResistancePlasticLaw &b, IceFlowLaw &i, EnthalpyConverter &e,
        const NCConfigVariable &c);

  SSAStrengthExtension *strength_extension;

  virtual ~SSA() { 
    PetscErrorCode ierr = deallocate();
    if (ierr != 0) {
      PetscPrintf(grid.com, "FATAL ERROR: SSAFD de-allocation failed.\n");
      PISMEnd();
    }
    delete strength_extension;
  }

  virtual PetscErrorCode init(PISMVars &vars);

  virtual PetscErrorCode update(bool fast);

  virtual PetscErrorCode set_initial_guess(IceModelVec2V &guess);

  virtual PetscErrorCode stdout_report(string &result);

  virtual void add_vars_to_output(string keyword, set<string> &result);
  virtual PetscErrorCode define_variables(set<string> vars, const NCTool &nc, nc_type nctype);
  virtual PetscErrorCode write_variables(set<string> vars, string filename);

  virtual void get_diagnostics(map<string, PISMDiagnostic*> &dict);

  using ShallowStressBalance::compute_principal_strain_rates;
  virtual PetscErrorCode compute_principal_strain_rates(
                IceModelVec2S &result_e1, IceModelVec2S &result_e2);

protected:
  virtual PetscErrorCode allocate();

  virtual PetscErrorCode deallocate();

  virtual PetscErrorCode solve()  = 0;

  virtual PetscErrorCode compute_driving_stress(IceModelVec2V &taud);

  virtual PetscErrorCode compute_basal_frictional_heating(IceModelVec2S &result);

  virtual PetscErrorCode compute_D2(IceModelVec2S &result);

  virtual PetscErrorCode compute_maximum_velocity();

  IceModelVec2Int *mask;
  IceModelVec2S *thickness, *tauc, *surface, *bed;
  IceModelVec2V taud, velocity_old;
  IceModelVec3 *enthalpy;

  string stdout_ssa;

  // objects used by the SSA solver (internally)
  DA  SSADA;                    // dof=2 DA (grid.da2 has dof=1)
  Vec SSAX;  // global vector for solution

  // profiling
  int event_ssa;
};

class SSA_taud : public PISMDiag<SSA>
{
public:
  SSA_taud(SSA *m, IceGrid &g, PISMVars &my_vars);
  PetscErrorCode compute(IceModelVec* &result);
};

#endif /* _SSA_H_ */