From b3816e9331835f92bbbfe102a7fdb9483fa04d6b Mon Sep 17 00:00:00 2001 From: Constantine Khrulev Date: Tue, 9 Jun 2020 08:56:46 -0800 Subject: [PATCH] Remove the old BP solver implementation --- .../blatter/BlatterStressBalance.cc | 484 -------- .../blatter/BlatterStressBalance.hh | 141 --- .../blatter/Blatter_implementation.c | 1078 ----------------- .../blatter/Blatter_implementation.h | 137 --- src/stressbalance/blatter/CMakeLists.txt | 2 - src/stressbalance/blatter/FE3DTools.h | 314 ----- 6 files changed, 2156 deletions(-) delete mode 100644 src/stressbalance/blatter/BlatterStressBalance.cc delete mode 100644 src/stressbalance/blatter/BlatterStressBalance.hh delete mode 100644 src/stressbalance/blatter/Blatter_implementation.c delete mode 100644 src/stressbalance/blatter/Blatter_implementation.h delete mode 100644 src/stressbalance/blatter/FE3DTools.h diff --git a/src/stressbalance/blatter/BlatterStressBalance.cc b/src/stressbalance/blatter/BlatterStressBalance.cc deleted file mode 100644 index 255042cc20..0000000000 --- a/src/stressbalance/blatter/BlatterStressBalance.cc +++ /dev/null @@ -1,484 +0,0 @@ -// Copyright (C) 2010-2016, 2019, 2020 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 3 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 - -#include "BlatterStressBalance.hh" -#include "pism/coupler/OceanModel.hh" -#include "pism/util/IceGrid.hh" -#include "pism/util/Vars.hh" -#include "pism/basalstrength/basal_resistance.hh" -#include "FE3DTools.h" -#include "pism/util/EnthalpyConverter.hh" -#include "pism/rheology/FlowLaw.hh" -#include "pism/rheology/FlowLawFactory.hh" -#include "pism/util/Vars.hh" -#include "pism/util/error_handling.hh" -#include "pism/util/pism_utilities.hh" -#include "pism/geometry/Geometry.hh" -#include "pism/util/node_types.hh" -#include "pism/util/Context.hh" - -namespace pism { -namespace stressbalance { - -/*! - * FIXMEs: - * - * We need to allow spatially-variable (and depth-dependent) ice hardness. - * - * We need to compute the volumetric strain heating. - * - */ - -//! C-wrapper for PISM's IceFlowLaw::viscosity(). -extern "C" -PetscErrorCode viscosity(void *ctx, double hardness, double gamma, - double *eta, double *deta) { - BlatterQ1Ctx *blatter_ctx = (BlatterQ1Ctx*)ctx; - BlatterStressBalance *blatter_stress_balance = (BlatterStressBalance*)blatter_ctx->extra; - - try { - blatter_stress_balance->flow_law()->effective_viscosity(hardness, gamma, eta, deta); - } catch (...) { - MPI_Comm com = blatter_stress_balance->grid()->ctx()->com(); - handle_fatal_errors(com); - SETERRQ(com, 1, "A PISM callback failed"); - } - return 0; -} - -//! C-wrapper for PISM's IceBasalResistancePlasticLaw::dragWithDerivative(). -extern "C" -PetscErrorCode drag(void *ctx, double tauc, double u, double v, - double *taud, double *dtaub) { - BlatterQ1Ctx *blatter_ctx = (BlatterQ1Ctx*)ctx; - BlatterStressBalance *blatter_stress_balance = (BlatterStressBalance*)blatter_ctx->extra; - - try { - blatter_stress_balance->sliding_law()->drag_with_derivative(tauc, u, v, taud, dtaub); - } catch (...) { - MPI_Comm com = blatter_stress_balance->grid()->ctx()->com(); - handle_fatal_errors(com); - SETERRQ(com, 1, "A PISM callback failed"); - } - return 0; -} - -/*! @brief A no-cost wrapper around 2D arrays. Hides the indexing order. */ -template -class PointerWrapper2D { - T **m_data; -public: - PointerWrapper2D() { - m_data = NULL; - } - - T*** address() { - return &m_data; - } - - T& operator()(int i, int j) { - // NOTE: this indexing order is important! - return m_data[i][j]; - } -}; - -/*! @brief A no-cost wrapper around 3D arrays. Hides the indexing order. */ -template -class PointerWrapper3D { - T ***m_data; -public: - PointerWrapper3D() { - m_data = NULL; - } - - T**** address() { - return &m_data; - } - - T& operator()(int i, int j, int k) { - // NOTE: this indexing order is important! - return m_data[i][j][k]; - } -}; - -BlatterStressBalance::BlatterStressBalance(IceGrid::ConstPtr grid, - EnthalpyConverter::Ptr e) - : ShallowStressBalance(grid), - m_ice_bottom_surface(grid, "ice_bottom_surface", WITHOUT_GHOSTS), - m_node_type(grid, "node_type", WITHOUT_GHOSTS), - m_u(grid, "uvel", WITH_GHOSTS), - m_v(grid, "vvel", WITH_GHOSTS), - m_strain_heating(grid, "strainheat", WITHOUT_GHOSTS), // never diff'ed in hor dirs - m_min_thickness(10.0) { - - Config::ConstPtr config = grid->ctx()->config(); - - int blatter_Mz = (int)config->get_number("stress_balance.blatter.Mz"); - m_da2 = grid->get_dm(1, (int)config->get_number("grid.max_stencil_width")); - - m_ctx.Lx = 2.0 * grid->Lx(); - m_ctx.Ly = 2.0 * grid->Ly(); - m_ctx.dirichlet_scale = 1.0; - m_ctx.rhog = config->get_number("constants.ice.density") * config->get_number("constants.standard_gravity"); - m_ctx.no_slip = PETSC_TRUE; // FIXME (at least make configurable) - m_ctx.nonlinear.viscosity = viscosity; - m_ctx.nonlinear.drag = drag; - m_ctx.extra = this; - initialize_Q12D(m_ctx.Q12D.chi, m_ctx.Q12D.dchi); - initialize_Q13D(m_ctx.Q13D.chi, m_ctx.Q13D.dchi); - - PetscErrorCode ierr = BlatterQ1_create(grid->com, *m_da2, blatter_Mz, - &this->m_ctx, this->m_snes.rawptr()); - PISM_CHK(ierr, "BlatterQ1_create"); - - m_u.set_attrs("diagnostic", "horizontal velocity of ice in the X direction", - "m s-1", "m year-1", "land_ice_x_velocity", 0); - - m_v.set_attrs("diagnostic", "horizontal velocity of ice in the Y direction", - "m s-1", "m year-1", "land_ice_y_velocity", 0); - - m_strain_heating.set_attrs("internal", - "rate of strain heating in ice (dissipation heating)", - "W m-3", "mW m-3", "", 0); - - - // Storage for u and v on the sigma vertical grid (for restarting) - // - // These IceModelVec3Custom instances use sigma levels computed here and so cannot be - // allocated in the member initializer list above - { - std::vector sigma(blatter_Mz); - double dz = 1.0 / (blatter_Mz - 1); - for (int i = 0; i < blatter_Mz; ++i) { - sigma[i] = i * dz; - } - sigma.back() = 1.0; - - std::map z_attrs = - {{"axis", "Z"}, - {"long_name", "scaled Z-coordinate in the ice (z_base=0, z_surface=1)"}, - {"units", "1"}, - {"positive", "up"}}; - - m_u_sigma.reset(new IceModelVec3Custom(grid, "uvel_sigma", "z_sigma", sigma, z_attrs)); - m_u_sigma->set_attrs("diagnostic", - "horizontal velocity of ice in the X direction on the sigma vertical grid", - "m s-1", "m year-1", "", 0); - - m_v_sigma.reset(new IceModelVec3Custom(grid, "vvel_sigma", "z_sigma", sigma, z_attrs)); - m_v_sigma->set_attrs("diagnostic", - "horizontal velocity of ice in the Y direction on the sigma vertical grid", - "m s-1", "m year-1", "", 0); - } - - { - rheology::FlowLawFactory ice_factory("stress_balance.blatter.", config, e); - ice_factory.remove(ICE_GOLDSBY_KOHLSTEDT); - - ice_factory.set_default(config->get_string("stress_balance.blatter.flow_law")); - - m_flow_law = ice_factory.create(); - } -} - -BlatterStressBalance::~BlatterStressBalance() { -} - -void BlatterStressBalance::init_impl() { - // empty -} - -void BlatterStressBalance::update(const Inputs &inputs, bool full_update) { - - assert(full_update); - - // setup - setup(inputs); - - // solve - PetscErrorCode ierr = SNESSolve(this->m_snes, PETSC_NULL, PETSC_NULL); - PISM_CHK(ierr, "SNESSolve"); - - // Transfer solution from the FEM mesh to the regular grid used in the rest - // of PISM and compute the vertically-averaged velocity. - transfer_velocity(inputs.geometry->ice_thickness); - - // Copy solution from the SNES to m_u_sigma and m_v_sigma. - copy_velocity(FROM_SNES_STORAGE); - - compute_volumetric_strain_heating(); -} - -/*! \brief Set up model parameters on the fine grid. */ -/*! - * This method expects bed_elevation, ice_thickness, and tauc to have width=1 ghosts. - * - * We should also compute ice hardness on the "sigma" grid here. - */ -void BlatterStressBalance::setup(const Inputs &inputs) { - PetscErrorCode ierr; - Vec param_vec; - PointerWrapper2D parameters; - DM da; - - initialize_ice_hardness(*inputs.enthalpy, inputs.geometry->ice_thickness); - - ierr = SNESGetDM(this->m_snes, &da); PISM_CHK(ierr, "SNESGetDM"); - - ierr = BlatterQ1_begin_2D_parameter_access(da, PETSC_FALSE, ¶m_vec, parameters.address()); - PISM_CHK(ierr, "BlatterQ1_begin_2D_parameter_access"); - - const IceModelVec2S - &bed = inputs.geometry->bed_elevation, - &thickness = inputs.geometry->ice_thickness, - &tauc = *inputs.basal_yield_stress; - - ice_bottom_surface(*inputs.geometry, m_ice_bottom_surface); - - compute_node_types(inputs.geometry->ice_thickness, m_min_thickness, m_node_type); - - IceModelVec::AccessList list{&bed, &thickness, &tauc, &m_ice_bottom_surface, - &m_node_type}; - - for (Points p(*m_grid); p; p.next()) { - const int i = p.i(), j = p.j(); - - parameters(i, j).ice_bottom = m_ice_bottom_surface(i, j); - - parameters(i, j).thickness = thickness(i,j); - - // fudge ice thickness (FIXME!!!) - if (thickness(i,j) < m_min_thickness) - parameters(i, j).thickness += m_min_thickness; - - parameters(i, j).tauc = tauc(i,j); - - parameters(i, j).node_type = m_node_type(i, j); - } - - ierr = BlatterQ1_end_2D_parameter_access(da, PETSC_FALSE, ¶m_vec, parameters.address()); - PISM_CHK(ierr, "BlatterQ1_end_2D_parameter_access"); -} - -//! Initialize ice hardness on the "sigma" grid. -void BlatterStressBalance::initialize_ice_hardness(const IceModelVec3 &enthalpy, - const IceModelVec2S &ice_thickness) { - PetscErrorCode ierr; - - PointerWrapper3D hardness; - unsigned int Mz_fem = static_cast(m_config->get_number("stress_balance.blatter.Mz")); - DM da; - Vec hardness_vec; - - ierr = SNESGetDM(this->m_snes, &da); PISM_CHK(ierr, "SNESGetDM"); - - ierr = BlatterQ1_begin_hardness_access(da, PETSC_FALSE, &hardness_vec, hardness.address()); - PISM_CHK(ierr, "BlatterQ1_begin_hardness_access"); - - const std::vector &zlevels = enthalpy.levels(); - - IceModelVec::AccessList list{&enthalpy, &ice_thickness}; - - for (Points p(*m_grid); p; p.next()) { - const int i = p.i(), j = p.j(); - - double thk = ice_thickness(i, j); - - // fudge ice thickness (FIXME!!!) - if (thk < m_min_thickness) - thk += m_min_thickness; - - double dz_fem = thk / (Mz_fem - 1); - const double *E = enthalpy.get_column(i, j); - - // compute ice hardness on the sigma grid - for (unsigned int k = 0; k < Mz_fem; ++k) { - double z_fem = k * dz_fem, - depth = thk - z_fem, - pressure = m_EC->pressure(depth), - E_local; - unsigned int k0 = m_grid->kBelowHeight(z_fem); - - const unsigned int Mz = m_grid->Mz(); - - if (k0 + 1 < Mz) { - double lambda = (z_fem - zlevels[k0]) / (zlevels[k0+1] - zlevels[k0]); - - E_local = (1.0 - lambda) * E[k0] + lambda * E[k0 + 1]; - } else { - // should never happen - E_local = E[Mz-1]; - } - - hardness(i, j, k) = m_flow_law->hardness(E_local, pressure); - } - } - - ierr = BlatterQ1_end_hardness_access(da, PETSC_FALSE, &hardness_vec, hardness.address()); - PISM_CHK(ierr, "BlatterQ1_end_hardness_access"); -} - -//! Transfer the velocity field from the FEM "sigma" grid onto PISM's grid. -/*! - * We also compute vertically-averaged ice velocity here. - */ -void BlatterStressBalance::transfer_velocity(const IceModelVec2S &ice_thickness) { - PetscErrorCode ierr; - - PointerWrapper3D U; - double *u_ij, *v_ij; - DM da; - Vec X; - unsigned int Mz_fem = static_cast(m_config->get_number("stress_balance.blatter.Mz")); - - ierr = SNESGetDM(this->m_snes, &da); PISM_CHK(ierr, "SNESGetDM"); - ierr = SNESGetSolution(this->m_snes, &X); PISM_CHK(ierr, "SNESGetSolution"); - - ierr = DMDAVecGetArray(da, X, U.address()); PISM_CHK(ierr, "DMDAVecGetArray"); - - IceModelVec::AccessList list{&m_u, &m_v, &ice_thickness, &m_velocity}; - - const unsigned int Mz = m_grid->Mz(); - const std::vector &zlevels = m_u.levels(); - - for (Points p(*m_grid); p; p.next()) { - const int i = p.i(), j = p.j(); - - u_ij = m_u.get_column(i, j); - v_ij = m_v.get_column(i, j); - - double thk = ice_thickness(i,j); - - // fudge ice thickness (FIXME!!!) - if (thk < m_min_thickness) - thk += m_min_thickness; - - double dz_fem = thk / (Mz_fem - 1); - - // compute vertically-averaged velocity using trapezoid rule - double ubar = 0, vbar = 0; - for (unsigned int k = 0; k < Mz_fem - 1; ++k) { - ubar += U(i, j, k).u + U(i, j, k+1).u; - vbar += U(i, j, k).v + U(i, j, k+1).v; - } - // finish the traperoidal rule (1/2 * dz) and compute the average: - m_velocity(i,j).u = ubar * (0.5*dz_fem) / thk; - m_velocity(i,j).v = vbar * (0.5*dz_fem) / thk; - - // compute 3D horizontal velocity - unsigned int current_level = 0; - for (unsigned int k = 0; k < Mz; ++k) { - - // find the FEM grid level just below the current PISM grid level - while ((current_level + 1) * dz_fem < zlevels[k]) - current_level++; - - if (current_level + 1 < Mz_fem) { - // use linear interpolation - double z0 = current_level * dz_fem, - lambda = (zlevels[k] - z0) / dz_fem; - - u_ij[k] = (U(i, j, current_level).u * (1 - lambda) + - U(i, j, current_level+1).u * lambda); - - v_ij[k] = (U(i, j, current_level).v * (1 - lambda) + - U(i, j, current_level+1).v * lambda); - - } else { - // extrapolate above the surface - u_ij[k] = U(i, j, Mz_fem-1).u; - v_ij[k] = U(i, j, Mz_fem-1).v; - } - - } // k-loop - } // loop over map-plane grid points - - ierr = DMDAVecRestoreArray(da, X, U.address()); PISM_CHK(ierr, "DMDAVecRestoreArray"); - - m_u.update_ghosts(); - m_v.update_ghosts(); -} - -//! Copy velocity from a dof=2 vector to special storage (to save it for re-starting). -void BlatterStressBalance::copy_velocity(Direction direction) { - PetscErrorCode ierr; - - PointerWrapper3D U; - DM da; - Vec X; - unsigned int Mz_fem = static_cast(m_config->get_number("stress_balance.blatter.Mz")); - - ierr = SNESGetDM(this->m_snes, &da); PISM_CHK(ierr, "SNESGetDM"); - ierr = SNESGetSolution(this->m_snes, &X); PISM_CHK(ierr, "SNESGetSolution"); - - ierr = DMDAVecGetArray(da, X, U.address()); PISM_CHK(ierr, "DMDAVecGetArray"); - - IceModelVec::AccessList list{m_u_sigma.get(), m_v_sigma.get()}; - - for (Points p(*m_grid); p; p.next()) { - const int i = p.i(), j = p.j(); - - double - *u_ij = m_u_sigma->get_column(i, j), - *v_ij = m_v_sigma->get_column(i, j); - - if (direction == FROM_SNES_STORAGE) { - for (unsigned int k = 0; k < Mz_fem; ++k) { - u_ij[k] = U(i, j, k).u; - v_ij[k] = U(i, j, k).v; - } - } else { - for (unsigned int k = 0; k < Mz_fem; ++k) { - U(i, j, k).u = u_ij[k]; - U(i, j, k).v = v_ij[k]; - } - } - } // loop over map-plane grid points - - ierr = DMDAVecRestoreArray(da, X, U.address()); PISM_CHK(ierr, "DMDAVecRestoreArray"); -} - -const IceModelVec3& BlatterStressBalance::velocity_u() const { - return m_u; -} - -const IceModelVec3& BlatterStressBalance::velocity_v() const { - return m_v; -} - -const IceModelVec3Custom& BlatterStressBalance::velocity_u_sigma() const { - return *m_u_sigma; -} - -const IceModelVec3Custom& BlatterStressBalance::velocity_v_sigma() const { - return *m_v_sigma; -} - -void BlatterStressBalance::compute_volumetric_strain_heating() { - // FIXME: implement this - m_strain_heating.set(0.0); -} - -//! \brief Produce a report string for the standard output. -std::string BlatterStressBalance::stdout_report() const { - return ""; -} - - -} // end of namespace stressbalance -} // end of namespace pism diff --git a/src/stressbalance/blatter/BlatterStressBalance.hh b/src/stressbalance/blatter/BlatterStressBalance.hh deleted file mode 100644 index 1da5b1c074..0000000000 --- a/src/stressbalance/blatter/BlatterStressBalance.hh +++ /dev/null @@ -1,141 +0,0 @@ -// Copyright (C) 2010-2016, 2018, 2019, 2020 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 3 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 _BLATTERSTRESSBALANCE_H_ -#define _BLATTERSTRESSBALANCE_H_ - -#include "pism/stressbalance/ShallowStressBalance.hh" -#include "pism/stressbalance/StressBalance.hh" // Inputs -#include "pism/util/IceGrid.hh" -#include "pism/util/iceModelVec3Custom.hh" - -#include "pism/util/petscwrappers/DM.hh" -#include "pism/util/petscwrappers/SNES.hh" - -#include "Blatter_implementation.h" - -namespace pism { -namespace stressbalance { - -//! \brief Blatter-Pattyn stress balance based on Jed Brown's PETSc -//! tutorial ex48.c (Brown et al. 2011). -/*! -Toy hydrostatic ice flow with multigrid in 3D - -Solves the hydrostatic (aka Blatter/Pattyn/First Order) equations for ice sheet -flow using multigrid. The ice uses a power-law rheology with Glen exponent 3 -(corresponds to p=4/3 in a p-Laplacian). - -The equations for horizontal velocity @f$ (u,v) @f$ are -@f{align*} - - [\eta (4 u_x + 2 v_y)]_x - [\eta (u_y + v_x)]_y - [\eta u_z]_z + \rho g s_x &= 0 \\ - - [\eta (4 v_y + 2 u_x)]_y - [\eta (u_y + v_x)]_x - [\eta v_z]_z + \rho g s_y &= 0 -@f} -where - @f[\eta = B/2 (\epsilon + \gamma)^{(p-2)/2}.@f] -is the nonlinear effective viscosity with regularization epsilon and hardness -parameter B, written in terms of the second invariant - -@f[\gamma = u_x^2 + v_y^2 + u_x v_y + \frac{1}{4} (u_y + v_x)^2 + \frac{1}{4} u_z^2 + \frac{1}{4} v_z^2.@f] - -The surface boundary conditions are the natural conditions, -corresponding to the "zero stress" condition (see [\ref RappazReist05]). The -basal boundary conditions are either no-slip, or a pseudo-plastic -sliding law (see IceBasalResistancePlasticLaw). - -In the code, the equations for @f$ (u,v) @f$ are multiplied through by @f$ 1/(\rho g) @f$ -so that residuals are @f$ O(1) @f$. - -The discretization is @f$ Q_1 @f$ finite elements, managed by a DA. The grid is never -distorted in the map @f$ (x,y) @f$ plane, but the bed and surface may be bumpy. -This is handled as usual in FEM, through the Jacobian of the coordinate -transformation from a reference element to the physical element. - -Since ice-flow is tightly coupled in the z-direction (within columns), the DA is -managed specially so that columns are never distributed, and are always -contiguous in memory. This amounts to reversing the meaning of X,Y,Z compared -to the DA's internal interpretation, and then indexing as `vec[i][j][k]`. The -exotic coarse spaces require 2D DAs which are made to use compatible domain -decomposition relative to the 3D DAs. - -Note that this implementation introduces two extra simplifications -compatible with the small bed slope assumption: - -- the code evaluating the integral corresponding to the basal boundary - condition assumes that the Jacobian of the map from the 2D reference - element is @f$ J = \frac{1}{4} \Delta x \times \Delta y @f$, which is - correct only if the bed is a horizontal plane. - -- it assumes that the horizontal ice velocity at the base approximates - the tangential basal ice velocity, which is also correct if the base - of the ice is horizontal. - -See the source code `$PETSC_DIR/src/snes/examples/tutorials/ex48.c` for -the original implementation. - */ -class BlatterStressBalance : public ShallowStressBalance { -public: - BlatterStressBalance(IceGrid::ConstPtr g, EnthalpyConverter::Ptr e); - virtual ~BlatterStressBalance(); - - virtual std::string stdout_report() const; - - virtual void update(const Inputs &inputs, bool full_update); - - const IceModelVec3& velocity_u() const; - const IceModelVec3& velocity_v() const; - - const IceModelVec3Custom& velocity_u_sigma() const; - const IceModelVec3Custom& velocity_v_sigma() const; - -protected: - void init_impl(); - - void transfer_velocity(const IceModelVec2S &ice_thickness); - void initialize_ice_hardness(const IceModelVec3 &enthalpy, - const IceModelVec2S &ice_thickness); - void setup(const Inputs &inputs); - - void compute_volumetric_strain_heating(); - - enum Direction {FROM_SNES_STORAGE, TO_SNES_STORAGE}; - - void copy_velocity(Direction); - - IceModelVec2S m_ice_bottom_surface; - IceModelVec2Int m_node_type; - - IceModelVec3 m_u, m_v, m_strain_heating; - - // u and v components on the "sigma" vertical grid - IceModelVec3Custom::Ptr m_u_sigma, m_v_sigma; - - BlatterQ1Ctx m_ctx; - petsc::SNES m_snes; - - petsc::DM::Ptr m_da2; - - double m_min_thickness; // FIXME: this should be used to set boundary conditions at ice margins - std::string m_stdout_blatter; -}; - -} // end of namespace stressbalance -} // end of namespace pism - -#endif /* _BLATTERSTRESSBALANCE_H_ */ - diff --git a/src/stressbalance/blatter/Blatter_implementation.c b/src/stressbalance/blatter/Blatter_implementation.c deleted file mode 100644 index 20fa612c23..0000000000 --- a/src/stressbalance/blatter/Blatter_implementation.c +++ /dev/null @@ -1,1078 +0,0 @@ -/* Copyright (C) 2013, 2014, 2015, 2016, 2019, 2020 Jed Brown and the PISM Authors - - 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 3 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 -*/ - -#include -#include -#include -#include -#include - -#include "FE3DTools.h" - -#include "Blatter_implementation.h" - -/*! \brief Compute the square of a number */ -static PetscScalar Sqr(PetscScalar a) {return a*a;} - -/*! \file Blatter_implementation.c - - This file contains the implementation of the \f$Q_1\f$ 3D FEM solver for - the Blatter-Pattyn stress balance system. - - Define the Blatter effective strain rate tensor \f$M\f$ - \f[ - M = - \left( - \begin{array}{lll} - 4 u_x + 2 v_y & u_y + v_x & u_z\\ - u_y + v_x & 2 u_x + 4 v_y & v_z - \end{array} - \right). - \f] - - Then the Blatter system of equations reads - \f[ - \nabla\cdot (\eta M) + \rho g \nabla s = 0, - \f] - or - \f{align*}{ - \left(\eta (4 u_x + 2 v_y)\right)_x + \left(\eta (u_y + v_x)\right)_y + \left(\eta u_z\right)_z + \rho g s_x &= 0,\\ - \left(\eta (2 u_x + 4 v_y)\right)_y + \left(\eta (u_y + v_x)\right)_x + \left(\eta v_z\right)_z + \rho g s_y &= 0, - \f} - where \f$\eta\f$ is the effective viscosity - \f[ - \eta = \frac B 2 \left(\gamma + \epsilon^2\right)^{\displaystyle\frac{1-n}{2n}} - \f] - and \f$B\f$ and \f$\epsilon\f$ are ice hardness and the regularizing parameter, respectively. - - Also, \f$\gamma\f$ is the second invariant of the strain rate tensor, defined by - \f[ - \gamma(u,v) = u_{x}^2 + v_{y}^2 + v_{y} \cdot u_{x} + \frac14{(u_{y}+v_{x})^2} + \frac14{u_{z}^2} + \frac14{v_{z}^2} - \f] - See compute_nonlinearity() for the computation of \f$\gamma\f$ and \f$\eta\f$. -*/ - -typedef PetscErrorCode (*DMDASNESJacobianLocal)(DMDALocalInfo*, void*, Mat, Mat, MatStructure*, void*); -typedef PetscErrorCode (*DMDASNESFunctionLocal)(DMDALocalInfo*, void*, void*, void*); - -static void compute_surface_gradient(PetscReal dchi[4][4][2], const PrmNode parameters[], - PetscReal dx, PetscReal dy, PetscReal ds[8][2]); - -static void compute_nodal_z_coordinates(const PrmNode parameters[], PetscInt k, PetscInt zm, PetscReal zn[]); - -static PetscErrorCode BlatterQ1_restriction_hook(DM fine, - Mat mrestrict, Vec rscale, Mat inject, - DM coarse, void *ctx); - - -/*! \brief Set up the DM and allocate storage for model parameters on the current grid level. - */ -static PetscErrorCode BlatterQ1_setup_level(BlatterQ1Ctx *ctx, DM dm) -{ - PetscErrorCode ierr; - PetscInt refinelevel, coarsenlevel, level, Mx, My, Mz, mx, my, stencil_width; - DMDAStencilType stencil_type; - const PetscInt *lx, *ly; - - /* Storage (DMDA and Vec) for 2D parameters */ - DM da2prm; - Vec parameters; - /* Storage (DMDA and Vec) for 3D ice hardness */ - DM da3; - Vec hardness; - - MPI_Comm comm; - ierr = PetscObjectGetComm((PetscObject)dm, &comm); CHKERRQ(ierr); - - PetscFunctionBegin; - ierr = DMDAGetInfo(dm, NULL, /* dimensions */ - &Mz, &My, &Mx, /* grid size */ - NULL, &my, &mx, /* number of processors in each direction */ - NULL, /* number of degrees of freedom */ - &stencil_width, - NULL, NULL, NULL, /* types of ghost nodes at the boundary */ - &stencil_type); CHKERRQ(ierr); - - ierr = DMDAGetOwnershipRanges(dm, NULL, &ly, &lx); CHKERRQ(ierr); - - ierr = DMGetRefineLevel(dm, &refinelevel); CHKERRQ(ierr); - ierr = DMGetCoarsenLevel(dm, &coarsenlevel); CHKERRQ(ierr); - level = refinelevel - coarsenlevel; - - /* number of parameters per map-plane location */ - int dof = sizeof(PrmNode)/sizeof(PetscScalar); - - ierr = DMDACreate2d(comm, - DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, - stencil_type, - My, Mx, my, mx, /* grid size */ - dof, stencil_width, - ly, /* number of nodes per processor*/ - lx, /* ditto */ - &da2prm); CHKERRQ(ierr); - -#if PETSC_VERSION_GE(3,8,0) - ierr = DMSetUp(da2prm); CHKERRQ(ierr); -#endif - - { - ierr = PetscPrintf(comm, - "Level %D domain size (m) %8.2g x %8.2g, num elements %3d x %3d x %3d (%8d), size (m) %g x %g\n", - level, ctx->Lx, ctx->Ly, Mx, My, Mz, Mx*My*Mz, ctx->Lx / Mx, ctx->Ly / My); CHKERRQ(ierr); - } - - ierr = DMCreateGlobalVector(da2prm, ¶meters); CHKERRQ(ierr); - - ierr = PetscObjectCompose((PetscObject)dm, "DMDA_2D", (PetscObject)da2prm); CHKERRQ(ierr); - ierr = PetscObjectCompose((PetscObject)dm, "DMDA_2D_Vec", (PetscObject)parameters); CHKERRQ(ierr); - - ierr = DMDestroy(&da2prm); CHKERRQ(ierr); - ierr = VecDestroy(¶meters); CHKERRQ(ierr); - - - ierr = DMDACreate3d(comm, - DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, - DMDA_STENCIL_BOX, - Mz, My, Mx, - 1, my, mx, /* number of processors in z, y, x directions. (Always *one* in the z-direction.) */ - 1, /* number of degrees of freedom per node (one) */ - 1, /* stencil width */ - NULL, ly, lx, /* number of nodes per processor */ - &da3); CHKERRQ(ierr); - -#if PETSC_VERSION_GE(3,8,0) - ierr = DMSetUp(da3); CHKERRQ(ierr); -#endif - - ierr = DMCreateGlobalVector(da3, &hardness); CHKERRQ(ierr); - - ierr = PetscObjectCompose((PetscObject)dm, "DMDA_3D", (PetscObject)da3); CHKERRQ(ierr); - ierr = PetscObjectCompose((PetscObject)dm, "DMDA_3D_Vec", (PetscObject)hardness); CHKERRQ(ierr); - - ierr = DMDestroy(&da3); CHKERRQ(ierr); - ierr = VecDestroy(&hardness); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - - -/*! \brief Create the restriction matrix. - * - * The result of this call is attached to `dm_fine` under `mat_name`. - * - * \param[in] dm_fine DM corresponding to the fine grid - * \param[in] dm_coarse DM corresponding to the coarse grid - * \param[in] dm_name name of the DM ("DMDA_2D" or "DMDA_3D") - * \param[in] mat_name name to use when attaching the interpolation matrix to `dm_fine` - */ -static PetscErrorCode BlatterQ1_create_restriction(DM dm_fine, DM dm_coarse, - const char dm_name[], - const char mat_name[]) { - PetscErrorCode ierr; - DM da2_fine, da2_coarse; - Mat mat; - - /* 1. get the DM for parameters from the fine grid DM */ - ierr = PetscObjectQuery((PetscObject)dm_fine, dm_name, - (PetscObject*)&da2_fine); CHKERRQ(ierr); - if (!da2_fine) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No %s composed with given DMDA", dm_name); - - /* 2. get the DM for parameters from the coarse grid DM */ - ierr = PetscObjectQuery((PetscObject)dm_coarse, dm_name, - (PetscObject*)&da2_coarse); CHKERRQ(ierr); - if (!da2_coarse) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No %s composed with given DMDA", dm_name); - - /* call DMCreateInterpolation */ - ierr = DMCreateInterpolation(da2_coarse, da2_fine, - &mat, PETSC_NULL); CHKERRQ(ierr); - - /* attach to the fine grid DM */ - ierr = PetscObjectCompose((PetscObject)dm_fine, mat_name, - (PetscObject)mat); CHKERRQ(ierr); - ierr = MatDestroy(&mat); CHKERRQ(ierr); - - return 0; -} - -/*! \brief Grid coarsening hook. - * - * This hook is called *once* when SNES sets up the next coarse level. - * - * This hook does three things: - * - Set up the DM for the newly created coarse level. - * - Set up the matrix type on the coarsest level to allow using - * direct solvers for the coarse problem. - * - Set up the interpolation matrix that will be used by the - * restriction hook to set model parameters on the new coarse level. - * - * See BlatterQ1_restriction_hook(). - */ -static PetscErrorCode BlatterQ1_coarsening_hook(DM dm_fine, DM dm_coarse, void *ctx) -{ - PetscErrorCode ierr; - BlatterQ1Ctx *blatter_ctx = (BlatterQ1Ctx*)ctx; - PetscInt rlevel, clevel; - - PetscFunctionBegin; - ierr = BlatterQ1_setup_level(blatter_ctx, dm_coarse); CHKERRQ(ierr); - - ierr = DMGetRefineLevel(dm_coarse, &rlevel); CHKERRQ(ierr); - ierr = DMGetCoarsenLevel(dm_coarse, &clevel); CHKERRQ(ierr); - if (rlevel-clevel == 0) { - ierr = DMSetMatType(dm_coarse, MATAIJ); CHKERRQ(ierr); - } - - ierr = DMCoarsenHookAdd(dm_coarse, BlatterQ1_coarsening_hook, BlatterQ1_restriction_hook, - ctx); CHKERRQ(ierr); - - /* create 2D interpolation */ - ierr = BlatterQ1_create_restriction(dm_fine, dm_coarse, - "DMDA_2D", "DMDA_2D_Restriction"); CHKERRQ(ierr); - - /* create 3D interpolation */ - ierr = BlatterQ1_create_restriction(dm_fine, dm_coarse, - "DMDA_3D", "DMDA_3D_Restriction"); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -/*! \brief Restrict model parameters from the `fine` grid onto the `coarse` grid. - * - * This function uses the restriction matrix created by BlatterQ1_coarsening_hook(). - */ -static PetscErrorCode BlatterQ1_restrict(DM fine, DM coarse, - const char dm_name[], - const char mat_name[], - const char vec_name[]) -{ - PetscErrorCode ierr; - Vec X_fine, X_coarse; - DM da2_fine, da2_coarse; - Mat mat; - - PetscFunctionBegin; - - /* get the restriction matrix from the fine grid DM */ - ierr = PetscObjectQuery((PetscObject)fine, mat_name, - (PetscObject*)&mat); CHKERRQ(ierr); - if (!mat) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_Restriction composed with given DMDA"); - - /* get the DMDA from the fine grid DM */ - ierr = PetscObjectQuery((PetscObject)fine, dm_name, - (PetscObject*)&da2_fine); CHKERRQ(ierr); - if (!da2_fine) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_Vec composed with given DMDA"); - - /* get the storage vector from the fine grid DM */ - ierr = PetscObjectQuery((PetscObject)fine, vec_name, - (PetscObject*)&X_fine); CHKERRQ(ierr); - if (!X_fine) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_Vec composed with given DMDA"); - - /* get the DMDA from the coarse grid DM */ - ierr = PetscObjectQuery((PetscObject)coarse, dm_name, - (PetscObject*)&da2_coarse); CHKERRQ(ierr); - if (!da2_coarse) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_Vec composed with given DMDA"); - - /* get the storage vector from the coarse grid DM */ - ierr = PetscObjectQuery((PetscObject)coarse, vec_name, - (PetscObject*)&X_coarse); CHKERRQ(ierr); - if (!X_coarse) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_Vec composed with given DMDA"); - - ierr = MatRestrict(mat, X_fine, X_coarse); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - - -/*! \brief Restrict model parameters to the grid corresponding to the next coarse level. - - This function computes the product - \f[ X_{\text{coarse}} = R X_{\text{fine}}, \f] - - extracting necessary information from DM objects managing fine and coarse grids. - - This is called once per SNESSolve(). - */ -static PetscErrorCode BlatterQ1_restriction_hook(DM fine, - Mat mrestrict, Vec rscale, Mat inject, - DM coarse, void *ctx) -{ - PetscErrorCode ierr; - - /* Get rid of "unused argument" warnings: */ - (void) mrestrict; - (void) rscale; - (void) inject; - (void) ctx; - - ierr = BlatterQ1_restrict(fine, coarse, - "DMDA_2D", "DMDA_2D_Restriction", "DMDA_2D_Vec"); CHKERRQ(ierr); - - ierr = BlatterQ1_restrict(fine, coarse, - "DMDA_3D", "DMDA_3D_Restriction", "DMDA_3D_Vec"); CHKERRQ(ierr); - - return 0; -} - -/*! \brief Get the pointer to the 2D array storing models parameters. - - The input argument da is the 3D DM for the current level, but we - need the 2D DM managing 2D parameters, so we need to extract it - first. - - \param[in] da the DM managed by the SNES object - - \param[in] local if PETSC_TRUE, get a pointer to a local Vec (for - use in residual and Jacobian evaluation code); if - PETSC_FALSE, get a pointer to a global Vec (to set - parameter values) - - \param[out] X_out pointer to the Vec we're accessing (for passing to - BlatterQ1_end_2D_parameter_access) - - \param[out] prm pointer to the array -*/ -PetscErrorCode BlatterQ1_begin_2D_parameter_access(DM da, PetscBool local, Vec *X_out, - PrmNode ***prm) -{ - PetscErrorCode ierr; - DM da2prm; - Vec X; - PetscFunctionBegin; - - ierr = PetscObjectQuery((PetscObject)da, "DMDA_2D", - (PetscObject*)&da2prm); CHKERRQ(ierr); - if (!da2prm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_2D composed with given DMDA"); - - ierr = PetscObjectQuery((PetscObject)da, "DMDA_2D_Vec", - (PetscObject*)&X); CHKERRQ(ierr); - if (!X) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_2D_Vec composed with given DMDA"); - - if (local == PETSC_TRUE) { - ierr = DMGetLocalVector(da2prm, X_out); CHKERRQ(ierr); - ierr = DMGlobalToLocalBegin(da2prm, X, INSERT_VALUES, *X_out); CHKERRQ(ierr); - ierr = DMGlobalToLocalEnd(da2prm, X, INSERT_VALUES, *X_out); CHKERRQ(ierr); - } else { - *X_out = X; - } - - ierr = DMDAVecGetArray(da2prm, *X_out, prm); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -/*! \brief Restore the 2D array of model parameters. - - See BlatterQ1_begin_2D_parameter_access() for details. -*/ -PetscErrorCode BlatterQ1_end_2D_parameter_access(DM da, PetscBool local, Vec *X_out, PrmNode ***prm) -{ - PetscErrorCode ierr; - DM da2prm; - PetscFunctionBegin; - - ierr = PetscObjectQuery((PetscObject)da, "DMDA_2D", - (PetscObject*)&da2prm); CHKERRQ(ierr); - if (!da2prm) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_2D composed with given DMDA"); - - ierr = DMDAVecRestoreArray(da2prm, *X_out, prm); CHKERRQ(ierr); - - if (local == PETSC_TRUE) { - ierr = DMRestoreLocalVector(da2prm, X_out); CHKERRQ(ierr); - } - - PetscFunctionReturn(0); -} - -/*! \brief Get the 3D array of ice hardness. - * - * This is similar to BlatterQ1_begin_2D_parameter_access(), but for ice - * hardness. - */ -PetscErrorCode BlatterQ1_begin_hardness_access(DM da, PetscBool local, Vec *X_out, PetscScalar ****hardness) -{ - PetscErrorCode ierr; - DM da3; - Vec X; - PetscFunctionBegin; - - ierr = PetscObjectQuery((PetscObject)da, "DMDA_3D", - (PetscObject*)&da3); CHKERRQ(ierr); - if (!da3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_3D composed with given DMDA"); - - ierr = PetscObjectQuery((PetscObject)da, "DMDA_3D_Vec", - (PetscObject*)&X); CHKERRQ(ierr); - if (!X) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_3D_Vec composed with given DMDA"); - - if (local == PETSC_TRUE) { - ierr = DMGetLocalVector(da3, X_out); CHKERRQ(ierr); - ierr = DMGlobalToLocalBegin(da3, X, INSERT_VALUES, *X_out); CHKERRQ(ierr); - ierr = DMGlobalToLocalEnd(da3, X, INSERT_VALUES, *X_out); CHKERRQ(ierr); - } else { - *X_out = X; - } - - ierr = DMDAVecGetArray(da3, *X_out, hardness); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -/*! \brief Restore the 3D array of ice hardness. - - See BlatterQ1_begin_hardness_access() for details. -*/ -PetscErrorCode BlatterQ1_end_hardness_access(DM da, PetscBool local, Vec *X_out, PetscScalar ****hardness) -{ - PetscErrorCode ierr; - DM da3; - PetscFunctionBegin; - - ierr = PetscObjectQuery((PetscObject)da, "DMDA_3D", - (PetscObject*)&da3); CHKERRQ(ierr); - if (!da3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, - "No DMDA_3D composed with given DMDA"); - - ierr = DMDAVecRestoreArray(da3, *X_out, hardness); CHKERRQ(ierr); - - if (local == PETSC_TRUE) { - ierr = DMRestoreLocalVector(da3, X_out); CHKERRQ(ierr); - } - - PetscFunctionReturn(0); -} - - -/*! \brief Set the initial guess. */ -/*! - * FIXME: we need a callback similar to drag() and viscosity(). - */ -static PetscErrorCode BlatterQ1_initial_guess(SNES snes, Vec X, void* ctx) -{ - PetscErrorCode ierr; - - /* Get rid of "unused argument" warnings: */ - (void) snes; - (void) ctx; - - PetscFunctionBegin; - - ierr = VecSet(X, 0.0); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -/*! \brief Compute the viscosity non-linearity and related quantities. - * - * This function is called once for each quadrature point. It uses - * nodal values of \f$u\f$, \f$v\f$ and the basis expansion of \f$u\f$ - * and \f$v\f$ to compute \f$u\f$, \f$v\f$ at the quadrature point - * - * \f[ u(q_j) = \sum_{i=1}^8\phi_i(q_j)\cdot u_j \f] - * - * \f[ v(q_j) = \sum_{i=1}^8\phi_i(q_j)\cdot v_j \f] - * - * Here \f$q_j\f$ is a quadrature point, \f$u_j\f$ is a nodal value of \f$u\f$. - * - * In addition to this, it uses the map from the reference element to a - * physical element to compute partial derivatives of u,v with respect to x,y. - * - * \c gamma is the second invariant \f$ \gamma \f$. - * - * \param[in] ctx Solver's "application context". Provides `ctx->nonlinear.viscosity()`. - * \param[in] velocity nodal values of horizontal velocity - * \param[in] phi values of global basis functions \f$\phi\f$ (at the current quadrature point). - * \param[in] dphi values of derivatives of global basis functions \f$\phi\f$ (with respect to \f$x,y,z\f$) - * \param[out] u,v components of the horizontal velocity (at the current quadrature point) - * \param[out] du,dv partial derivatives of `u` and `v` - * \param[out] eta effective viscosity \f$\eta\f$ (at the current quadrature point) - * \param[out] deta derivative of eta with respect to gamma \f$\frac{\partial \eta}{\partial \gamma}\f$ - */ -static PetscErrorCode compute_nonlinearity(BlatterQ1Ctx *ctx, - const Node velocity[restrict], - const PetscReal hardness[restrict], - const PetscReal phi[restrict], - PetscReal dphi[restrict][3], - PetscScalar *restrict u, - PetscScalar *restrict v, - PetscScalar du[restrict], - PetscScalar dv[restrict], - PetscReal *eta, - PetscReal *deta) { - PetscErrorCode ierr; - PetscInt p, q; - PetscScalar second_invariant, hardness_q; - - /* Compute values of velocity components and their derivatives at the current quadrature - point. */ - du[0] = du[1] = du[2] = 0; - dv[0] = dv[1] = dv[2] = 0; - *u = 0; - *v = 0; - hardness_q = 0.0; - for (p = 0; p < 8; p++) { - *u += phi[p] * velocity[p].u; - *v += phi[p] * velocity[p].v; - hardness_q += phi[p] * hardness[p]; - for (q = 0; q < 3; q++) { - du[q] += dphi[p][q] * velocity[p].u; - dv[q] += dphi[p][q] * velocity[p].v; - } - } - - second_invariant = (Sqr(du[0]) + Sqr(dv[1]) + du[0]*dv[1] + - 0.25*Sqr(du[1] + dv[0]) + - 0.25*Sqr(du[2]) + - 0.25*Sqr(dv[2])); - - ierr = ctx->nonlinear.viscosity(ctx, hardness_q, second_invariant, eta, deta); CHKERRQ(ierr); - return 0; -} - - -/*! \brief Evaluate the residual. - * - * - * FIXME: I need to document this. - */ -static PetscErrorCode BlatterQ1_residual_local(DMDALocalInfo *info, Node ***velocity, Node ***residual, - BlatterQ1Ctx *ctx) -{ - PetscInt xs, ys, xm, ym, zm, i, j, k, q, l; - PetscReal dx, dy; - PrmNode **prm; - Vec prm_local; - PetscErrorCode ierr; - Vec hardness_local; - PetscReal ***hardness; - - PetscFunctionBegin; - xs = info->zs; - ys = info->ys; - xm = info->zm; - ym = info->ym; - zm = info->xm; - dx = ctx->Lx / info->mz; /* grid spacing in the x direction */ - dy = ctx->Ly / info->my; /* grid spacing in the y direction */ - ierr = BlatterQ1_begin_2D_parameter_access(info->da, PETSC_TRUE, - &prm_local, &prm); CHKERRQ(ierr); - - ierr = BlatterQ1_begin_hardness_access(info->da, PETSC_TRUE, - &hardness_local, &hardness); CHKERRQ(ierr); - - for (i = xs; i < xs + xm; i++) { - for (j = ys; j < ys + ym; j++) { - PrmNode parameters[4]; /* 4 nodes (in 2D) */ - PetscReal ds[8][2]; /* 8 quadrature points, 2 dimensions */ - get_nodal_values_2d(prm, i, j, parameters); - - compute_surface_gradient(ctx->Q12D.dchi, parameters, dx, dy, ds); - - /* loop over elements in the column */ - for (k = 0; k < zm - 1; k++) { - PetscInt ls = 0; /* starting index */ - Node U_nodal[8], *R_nodal[8]; - PetscReal hardness_nodal[8]; - PetscReal zn[8], etabase = 0; - - /* Compute z-coordinates of element nodes: */ - compute_nodal_z_coordinates(parameters, k, zm, zn); - - /* Get nodal values of velocity components: */ - get_nodal_values_3d(velocity, i, j, k, U_nodal); - get_nodal_values_3d(hardness, i, j, k, hardness_nodal); - - /* Get pointers to residual components corresponding to element nodes: */ - get_pointers_to_nodal_values_3d(residual, i, j, k, R_nodal); - - if (ctx->no_slip && k == 0) { - for (l = 0; l < 4; l++) { - U_nodal[l].u = 0; - U_nodal[l].v = 0; - } - /* The first 4 basis functions lie on the bottom layer, so their contribution is - * exactly 0 and we can skip them */ - ls = 4; - } - - for (q = 0; q < 8; q++) { /* for all quadrature points... */ - PetscReal grad_z[3], phi[8], dphi[8][3], W, eta, deta; - PetscScalar du[3], dv[3], u, v; - - /* Compute various quantities at this quadrature point in the *physical* element */ - - compute_z_gradient(ctx->Q13D.dchi[q], zn, /* inputs (derivatives of shape functions, geometry) */ - grad_z); /* output (partial derivatives of z with respect to xi,eta,zeta) */ - - compute_element_info(ctx->Q13D.chi, ctx->Q13D.dchi, q, dx, dy, grad_z, /* inputs */ - phi, dphi, &W); /* outputs (values of shape functions, their derivatives, weights */ - - ierr = compute_nonlinearity(ctx, U_nodal, hardness_nodal, phi, dphi, /* inputs */ - &u, &v, du, dv, &eta, &deta); /* outputs u,v, partial derivatives of u,v, - effective viscosity (eta), derivative of eta with respect to gamma */ - CHKERRQ(ierr); - - W /= ctx->rhog; /* scales residuals to be O(1) */ - - if (q == 0) { - etabase = eta; - } - - for (l = ls; l < 8; l++) { /* test functions */ - const PetscReal *dp = dphi[l]; - R_nodal[l]->u += dp[0]*W*eta*(4.0*du[0] + 2.0*dv[1]) + dp[1]*W*eta*(du[1] + dv[0]) + dp[2]*W*eta*du[2] + phi[l]*W*ctx->rhog*ds[q][0]; - R_nodal[l]->v += dp[1]*W*eta*(2.0*du[0] + 4.0*dv[1]) + dp[0]*W*eta*(du[1] + dv[0]) + dp[2]*W*eta*dv[2] + phi[l]*W*ctx->rhog*ds[q][1]; - } - } /* q-loop */ - - if (k == 0) { /* we are on a bottom face */ - if (ctx->no_slip) { - /* Note: Non-Galerkin coarse grid operators are very sensitive to the scaling of Dirichlet boundary - * conditions. After shenanigans above, etabase contains the effective viscosity at the closest quadrature - * point to the bed. We want the diagonal entry in the Dirichlet condition to have similar magnitude to the - * diagonal entry corresponding to the adjacent node. The fundamental scaling of the viscous part is in - * diagu, diagv below. This scaling is easy to recognize by considering the finite difference operator after - * scaling by element size. The no-slip Dirichlet condition is scaled by this factor, and also in the - * assembled matrix (see the similar block in BlatterQ1_Jacobian_local). - * - * Note that the residual at this Dirichlet node is linear in the state at this node, but also depends - * (nonlinearly in general) on the neighboring interior nodes through the local viscosity. This will make - * a matrix-free Jacobian have extra entries in the corresponding row. We assemble only the diagonal part, - * so the solution will exactly satisfy the boundary condition after the first linear iteration. - */ - const PetscReal dz = PetscRealPart(parameters[0].thickness) / (zm - 1.0); - const PetscScalar - diagu = 2*etabase / ctx->rhog*(dx*dy / dz + dx*dz / dy + 4*dy*dz / dx), - diagv = 2*etabase / ctx->rhog*(dx*dy / dz + 4*dx*dz / dy + dy*dz / dx); - R_nodal[0]->u = ctx->dirichlet_scale*diagu*velocity[i][j][k].u; - R_nodal[0]->v = ctx->dirichlet_scale*diagv*velocity[i][j][k].v; - } else { /* Integrate over bottom face to apply boundary condition */ - - for (q = 0; q < 4; q++) { /* for each quadrature point on the basal face... */ - /* Note: W below can be thought of as an approximation - of the Jacobian*weight using the small bed slope - assumption. (It *is* correct in the flat bed case.) - */ - const PetscReal W = 0.25*dx*dy / ctx->rhog, - *phi = ctx->Q12D.chi[q]; - - PetscScalar u = 0, v = 0, tauc = 0; - PetscReal beta; /* basal drag coefficient; tau_{b,x} = beta*u; tau_{b,y} = beta*v */ - for (l = 0; l < 4; l++) { - u += phi[l]*U_nodal[l].u; - v += phi[l]*U_nodal[l].v; - tauc += phi[l]*parameters[l].tauc; - } - - ierr = ctx->nonlinear.drag(ctx, PetscRealPart(tauc), u, v, - &beta, NULL); CHKERRQ(ierr); - - for (l = 0; l < 4; l++) { - R_nodal[ls + l]->u += phi[l]*W*(beta*u); - R_nodal[ls + l]->v += phi[l]*W*(beta*v); - } - } /* end of the quadrature loop */ - - } /* end of the generic basal boundary condition */ - - } /* end of "if (k == 0)" */ - } /* k-loop */ - } /* j-loop */ - } /* i-loop */ - - ierr = BlatterQ1_end_hardness_access(info->da, PETSC_TRUE, - &hardness_local, &hardness); CHKERRQ(ierr); - - ierr = BlatterQ1_end_2D_parameter_access(info->da, PETSC_TRUE, &prm_local, &prm); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -/*! \brief Evaluate the Jacobian. - * - * FIXME: I need to document this. - */ -static PetscErrorCode BlatterQ1_Jacobian_local(DMDALocalInfo *info, Node ***velocity, Mat A, Mat B, - BlatterQ1Ctx *ctx) { - PetscInt xs, ys, xm, ym, zm, i, j, k, q, l, ll; - PetscReal dx, dy; - PrmNode **prm; - Vec prm_local; - PetscErrorCode ierr; - Vec hardness_local; - PetscReal ***hardness; - - (void)A; - - PetscFunctionBegin; - xs = info->zs; - ys = info->ys; - xm = info->zm; - ym = info->ym; - zm = info->xm; - dx = ctx->Lx / info->mz; - dy = ctx->Ly / info->my; - - ierr = MatZeroEntries(B); CHKERRQ(ierr); - - ierr = BlatterQ1_begin_2D_parameter_access(info->da, PETSC_TRUE, - &prm_local, &prm); CHKERRQ(ierr); - - ierr = BlatterQ1_begin_hardness_access(info->da, PETSC_TRUE, - &hardness_local, &hardness); CHKERRQ(ierr); - - for (i = xs; i < xs + xm; i++) { - for (j = ys; j < ys + ym; j++) { - PrmNode parameters[4]; /* 4 nodes */ - get_nodal_values_2d(prm, i, j, parameters); - - for (k = 0; k < zm - 1; k++) { - PetscInt ls = 0; - Node U_nodal[8]; - PetscReal hardness_nodal[8]; - PetscReal zn[8], etabase = 0; - PetscScalar Ke[8*2][8*2]; - - PetscMemzero(Ke, sizeof(Ke)); - - compute_nodal_z_coordinates(parameters, k, zm, zn); - - get_nodal_values_3d(velocity, i, j, k, U_nodal); - get_nodal_values_3d(hardness, i, j, k, hardness_nodal); - - /* Ensure that values of velocity components at Dirichlet - (no-slip) nodes are exactly equal to Dirichlet B.C. values. */ - if (ctx->no_slip && k == 0) { - for (l = 0; l < 4; l++) { - U_nodal[l].u = 0.0; - U_nodal[l].v = 0.0; - } - ls = 4; - } - - for (q = 0; q < 8; q++) { - PetscReal grad_z[3], phi[8], dphi[8][3], W, eta, deta; - PetscScalar du[3], dv[3], u, v; - - /* Compute the gradient of z */ - compute_z_gradient(ctx->Q13D.dchi[q], zn, /* inputs */ - grad_z); /* output */ - - /* compute values of shape functions, their derivatives, and - quadrature weights at a quadrature point q. */ - compute_element_info(ctx->Q13D.chi, ctx->Q13D.dchi, q, dx, dy, grad_z, /* inputs */ - phi, dphi, &W); /* outputs */ - - /* Compute u,v, their derivatives, plus effective viscosity and its - derivative with respect to gamma. */ - ierr = compute_nonlinearity(ctx, U_nodal, hardness_nodal, phi, dphi, /* inputs */ - &u, &v, du, dv, &eta, &deta); /* outputs */ - CHKERRQ(ierr); - - W /= ctx->rhog; /* residuals are scaled by this factor */ - - /* Store eta at the quadrature point 0 (i.e. at a location near the bottom face of - * the element). It is used to scale the diagonal entries at the base in the - * no-slip case. */ - if (q == 0) { - etabase = eta; - } - - for (l = ls; l < 8; l++) { /* test functions */ - const PetscReal *restrict dp = dphi[l]; - - for (ll = l; ll < 8; ll++) { - const PetscReal *restrict dpl = dphi[ll]; - - /* The analytic Jacobian in nice, easy-to-read form */ - { - PetscScalar dgdu, dgdv; - /* Compute derivatives of gamma with respect to u and v. */ - dgdu = 2.0*du[0]*dpl[0] + dv[1]*dpl[0] + 0.5*(du[1] + dv[0])*dpl[1] + 0.5*du[2]*dpl[2]; - dgdv = 2.0*dv[1]*dpl[1] + du[0]*dpl[1] + 0.5*(du[1] + dv[0])*dpl[0] + 0.5*dv[2]*dpl[2]; - /* Picard part */ - Ke[l*2 + 0][ll*2 + 0] += dp[0]*W*eta*4.0*dpl[0] + dp[1]*W*eta*dpl[1] + dp[2]*W*eta*dpl[2]; - Ke[l*2 + 0][ll*2 + 1] += dp[0]*W*eta*2.0*dpl[1] + dp[1]*W*eta*dpl[0]; - Ke[l*2 + 1][ll*2 + 0] += dp[1]*W*eta*2.0*dpl[0] + dp[0]*W*eta*dpl[1]; - Ke[l*2 + 1][ll*2 + 1] += dp[1]*W*eta*4.0*dpl[1] + dp[0]*W*eta*dpl[0] + dp[2]*W*eta*dpl[2]; - /* extra Newton terms */ - Ke[l*2 + 0][ll*2 + 0] += dp[0]*W*deta*dgdu*(4.0*du[0] + 2.0*dv[1]) + dp[1]*W*deta*dgdu*(du[1] + dv[0]) + dp[2]*W*deta*dgdu*du[2]; - Ke[l*2 + 0][ll*2 + 1] += dp[0]*W*deta*dgdv*(4.0*du[0] + 2.0*dv[1]) + dp[1]*W*deta*dgdv*(du[1] + dv[0]) + dp[2]*W*deta*dgdv*du[2]; - Ke[l*2 + 1][ll*2 + 0] += dp[1]*W*deta*dgdu*(4.0*dv[1] + 2.0*du[0]) + dp[0]*W*deta*dgdu*(du[1] + dv[0]) + dp[2]*W*deta*dgdu*dv[2]; - Ke[l*2 + 1][ll*2 + 1] += dp[1]*W*deta*dgdv*(4.0*dv[1] + 2.0*du[0]) + dp[0]*W*deta*dgdv*(du[1] + dv[0]) + dp[2]*W*deta*dgdv*dv[2]; - } - } /* ll-loop (trial functions) */ - } /* l-loop (test functions) */ - } /* q-loop */ - - /* on a bottom face */ - if (k == 0) { - if (ctx->no_slip) { - const PetscReal dz = PetscRealPart(parameters[0].thickness) / (zm - 1); - const PetscScalar diagu = 2*etabase / ctx->rhog*(dx*dy / dz + dx*dz / dy + 4*dy*dz / dx), - diagv = 2*etabase / ctx->rhog*(dx*dy / dz + 4*dx*dz / dy + dy*dz / dx); - Ke[0][0] = ctx->dirichlet_scale*diagu; - Ke[1][1] = ctx->dirichlet_scale*diagv; - } else { - - for (q = 0; q < 4; q++) { - /* Note: W below can be thought of as an approximation - of the Jacobian*weight using the small bed slope - assumption. (It *is* correct in the flat bed case.) - */ - const PetscReal W = 0.25*dx*dy / ctx->rhog, - *phi = ctx->Q12D.chi[q]; - PetscScalar u = 0, v = 0, tauc = 0; - PetscReal beta, /* basal drag coefficient; tau_{b,x} = beta*u; tau_{b,y} = beta*v */ - dbeta; /* derivative of beta with respect to alpha = 1/2 * (u*u + v*v) */ - - /* Compute u, v, \tau_c at a quadrature point on the bottom face using basis expansions: */ - for (l = 0; l < 4; l++) { - u += phi[l]*U_nodal[l].u; - v += phi[l]*U_nodal[l].v; - tauc += phi[l]*parameters[l].tauc; - } - - /* Compute the friction coefficient at this quadrature point: */ - ierr = ctx->nonlinear.drag(ctx, PetscRealPart(tauc), u, v, - &beta, &dbeta); CHKERRQ(ierr); - - for (l = 0; l < 4; l++) { - const PetscReal pp = phi[l]; - for (ll = 0; ll < 4; ll++) { - const PetscReal ppl = phi[ll]; - Ke[l*2 + 0][ll*2 + 0] += pp*W*beta*ppl + pp*W*dbeta*u*u*ppl; - Ke[l*2 + 0][ll*2 + 1] += pp*W*dbeta*u*v*ppl; - Ke[l*2 + 1][ll*2 + 0] += pp*W*dbeta*v*u*ppl; - Ke[l*2 + 1][ll*2 + 1] += pp*W*beta*ppl + pp*W*dbeta*v*v*ppl; - } /* ll-loop (trial functions) */ - } /* l-loop (test functions) */ - } /* q-loop */ - - } /* generic basal boundary condition */ - } /* end of "if (k == 0)" */ - - /* Set matrix values: */ - { - const MatStencil rc[8] = {{i, j, k, 0}, {i + 1, j, k, 0}, {i + 1, j + 1, k, 0}, {i, j + 1, k, 0}, - {i, j, k + 1, 0}, {i + 1, j, k + 1, 0}, {i + 1, j + 1, k + 1, 0}, {i, j + 1, k + 1, 0}}; - { - /* fill in lower-triangular part, this is really cheap compared to computing the entries */ - for (l = 0; l < 8; l++) { - for (ll = l + 1; ll < 8; ll++) { - Ke[ll*2 + 0][l*2 + 0] = Ke[l*2 + 0][ll*2 + 0]; - Ke[ll*2 + 1][l*2 + 0] = Ke[l*2 + 0][ll*2 + 1]; - Ke[ll*2 + 0][l*2 + 1] = Ke[l*2 + 1][ll*2 + 0]; - Ke[ll*2 + 1][l*2 + 1] = Ke[l*2 + 1][ll*2 + 1]; - } - } - ierr = MatSetValuesBlockedStencil(B, 8, rc, 8, rc, &Ke[0][0], ADD_VALUES); CHKERRQ(ierr); - } - } - - } /* k-loop */ - } /* j-loop */ - } /* i-loop */ - - ierr = BlatterQ1_end_hardness_access(info->da, PETSC_TRUE, - &hardness_local, &hardness); CHKERRQ(ierr); - - ierr = BlatterQ1_end_2D_parameter_access(info->da, PETSC_TRUE, &prm_local, &prm); CHKERRQ(ierr); - - ierr = MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); - ierr = MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); - ierr = MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE); CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -/*! \brief Allocate the fine grid 3D DMDA and set up the SNES object. - - \param[in] com communicator - \param[in] pism_da PISM-side 2D DMDA (used to get grid info) - \param[in] Mz number of vertical levels - \param[in] ctx the application context - \param[out] result SNES object that will be used with SNESSolve later -*/ -PetscErrorCode BlatterQ1_create(MPI_Comm com, DM pism_da, - PetscInt Mz, - BlatterQ1Ctx *ctx, SNES *result) { - PetscErrorCode ierr; - PetscInt dim, Mx, My, Nx, Ny; - DM da; - const PetscInt *lx, *ly; - - ierr = DMDAGetInfo(pism_da, - &dim, - &Mx, - &My, - NULL, /* Mz */ - &Nx, /* number of processors in y-direction */ - &Ny, /* number of processors in x-direction */ - NULL, /* ditto, z-direction */ - NULL, /* number of degrees of freedom per node */ - NULL, /* stencil width */ - NULL, NULL, NULL, /* types of ghost nodes at the boundary */ - NULL); CHKERRQ(ierr); /* stencil type */ - assert(dim == 2); - - ierr = DMDAGetOwnershipRanges(pism_da, &lx, &ly, NULL); - - ierr = DMDACreate3d(com, - DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, - DMDA_STENCIL_BOX, - Mz, My, Mx, - 1, Ny, Nx, /* number of processors in z, y, x directions. (Always *one* in the z-direction.) */ - sizeof(Node)/sizeof(PetscScalar), - 1, /* stencil width */ - NULL, ly, lx, /* number of nodes per processor */ - &da); CHKERRQ(ierr); - -#if PETSC_VERSION_GE(3,8,0) - ierr = DMSetUp(da); CHKERRQ(ierr); -#endif - - ierr = DMDASetFieldName(da, 0, "x-velocity"); CHKERRQ(ierr); - ierr = DMDASetFieldName(da, 1, "y-velocity"); CHKERRQ(ierr); - - ierr = BlatterQ1_setup_level(ctx, da); CHKERRQ(ierr); - - /* ADD_VALUES, because BlatterQ1_residual_local contributes to ghosted values. */ - ierr = DMDASNESSetFunctionLocal(da, ADD_VALUES, - (DMDASNESFunction)BlatterQ1_residual_local, - ctx); CHKERRQ(ierr); - ierr = DMDASNESSetJacobianLocal(da, - (DMDASNESJacobian)BlatterQ1_Jacobian_local, - ctx); CHKERRQ(ierr); - - ierr = DMCoarsenHookAdd(da, BlatterQ1_coarsening_hook, - BlatterQ1_restriction_hook, ctx); CHKERRQ(ierr); - /* FIXME: do we need a refinement hook? */ - - ierr = DMSetApplicationContext(da, ctx); CHKERRQ(ierr); - - ierr = SNESCreate(com, result); CHKERRQ(ierr); - ierr = SNESSetDM(*result, da); CHKERRQ(ierr); - ierr = DMDestroy(&da); CHKERRQ(ierr); - - ierr = SNESSetComputeInitialGuess(*result, BlatterQ1_initial_guess, ctx); CHKERRQ(ierr); - ierr = SNESSetFromOptions(*result); CHKERRQ(ierr); - - return 0; -} - -/*! \brief Compute the surface gradient at all 4 quadrature points (in the map plane). - - Note that there are 8 quadrature points in a 3D hexahedral element, but since - surface elevation does not depend on \f$z\f$, values at the first 4 points (bottom - "level") are all we need. - - Using the 2D \f$Q_1\f$ basis expansion for the function \f$s(x,y)\f$, - - \f{eqnarray}{ - \frac{\partial s}{\partial x}(q_j) &=& \sum_{i=1}^4 \frac{\partial \phi_i}{\partial x}(q_j) s_i,\\ - \frac{\partial s}{\partial y}(q_j) &=& \sum_{i=1}^4 \frac{\partial \phi_i}{\partial y}(q_j) s_i. - \f} - - Here \f$s_i\f$ are nodal values of the surface elevation, \f$\phi_i\f$ are 2D - \f$Q_1\f$ global basis functions and \f$q_j\f$ are quadrature points. - - Now, if the grid has equal spacing in \f$x\f$ and \f$y\f$ directions, then - - \f{eqnarray}{ - \frac{\partial \phi_i}{\partial x} &=& \frac{2}{\Delta x} \frac{\partial \chi_i}{\partial \xi}\\ - \frac{\partial \phi_i}{\partial y} &=& \frac{2}{\Delta y} \frac{\partial \chi_i}{\partial \eta} - \f} - - This gives - \f{eqnarray}{ - \frac{\partial s}{\partial x}(q_j) &=& \frac{2}{\Delta x} \sum_{i=1}^4 \frac{\partial \chi_i}{\partial x}(q_j) s_i,\\ - \frac{\partial s}{\partial y}(q_j) &=& \frac{2}{\Delta y} \sum_{i=1}^4 \frac{\partial \chi_i}{\partial y}(q_j) s_i. - \f} - - Note also that 3D \f$Q_1\f$ element basis functions evaluated at an element face - reduce to 2D element basis functions, and so this is equivalent to assuming - that \f$s(x,y,z) = s(x,y)\f$ for all \f$z\f$ and using 3D \f$Q_1\f$ basis expansion - in this computation. - - \param[in] dchi values of derivatives of 2D element basis functions \f$\chi\f$ at quadrature points - \param[in] parameters 2D parameters at element nodes - \param[in] dx,dy grid spacing in x and y directions - \param[out] ds values of the surface gradient -*/ -static void compute_surface_gradient(PetscReal dchi[4][4][2], - const PrmNode parameters[], PetscReal dx, PetscReal dy, - PetscReal ds[8][2]) { - PetscInt i, q; - - /* loop over quadrature points */ - for (q = 0; q < 4; ++q) { - ds[q][0] = 0.0; - ds[q][1] = 0.0; - - /* loop over basis functions */ - for (i = 0; i < 4; ++i) { - double surface_elevation = parameters[i].ice_bottom + parameters[i].thickness; - ds[q][0] += dchi[q][i][0] * surface_elevation; - ds[q][1] += dchi[q][i][1] * surface_elevation; - } - - /* convert derivatives with respect to zeta and eta into - derivatives with respect to x and y */ - ds[q][0] *= 2.0/dx; - ds[q][1] *= 2.0/dy; - - /* In the 2x2x2 tensor-product quadrature the other four - quadrature points get the same values of ds/dx and ds/dy - (because s does not depend on z). */ - ds[q + 4][0] = ds[q][0]; - ds[q + 4][1] = ds[1][1]; - } -} - -/*! \brief Compute z-coordinates of the nodes of a hexahedral element. - - The bottom (top) mesh surface follows the bottom (top) surface of the ice. - Each "column" of nodes is equally-spaced, independently from others. - - \param[in] parameters 2D parameters at element nodes - \param[in] k index of the vertical level - \param[in] zm total number of vertical levels - \param[out] zn z-coordinates of element nodes -*/ -static void compute_nodal_z_coordinates(const PrmNode parameters[], PetscInt k, PetscInt zm, PetscReal zn[]) { - const double zm1 = zm - 1; - const double kk = k; - - /* loop over quadrature points */ - for (int i = 0; i < 4; i++) { - zn[i] = parameters[i].ice_bottom + parameters[i].thickness * kk / zm1; - zn[i+4] = parameters[i].ice_bottom + parameters[i].thickness * (kk + 1) / zm1; - } -} diff --git a/src/stressbalance/blatter/Blatter_implementation.h b/src/stressbalance/blatter/Blatter_implementation.h deleted file mode 100644 index afe29afdc1..0000000000 --- a/src/stressbalance/blatter/Blatter_implementation.h +++ /dev/null @@ -1,137 +0,0 @@ -/* Copyright (C) 2013, 2015, 2016, 2020 Jed Brown and the PISM Authors - * - * 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 3 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 _BLATTER_IMPLEMENTATION_H_ -#define _BLATTER_IMPLEMENTATION_H_ - -#include - -/*! This file contains declarations that need to be visible to the BlatterStressBalance C++ class. - * - */ - -#ifdef __cplusplus -extern "C" { -#endif - -/*! Storage for components of the horizontal ice velocity (at each node of the - 3D grid).*/ -typedef struct { - PetscScalar u, v; -} Node; - -/*! Storage for model parameters (at each node of the 2D grid). */ -typedef struct { - /*! elevation of the bottom surface of the ice */ - PetscScalar ice_bottom; - /*! thickness */ - PetscScalar thickness; - /*! till yield stress */ - PetscScalar tauc; - /*! node type */ - PetscScalar node_type; -} PrmNode; - -typedef struct { - /*! full width of the model domain in the x-direction */ - PetscReal Lx; - /*! full width of the model domain in the y-direction */ - PetscReal Ly; - /*! scaling for Dirichlet boundary conditions */ - PetscReal dirichlet_scale; - /*! ice density times standard gravity */ - PetscReal rhog; - /*! whether to use the no slip condition at the base */ - PetscBool no_slip; - struct { - /*! Function evaluating effective viscosity as a function of ice - hardness and the second invariant \f$ \gamma \f$ - - \param[in] ctx pointer to this context (BlatterQ1Ctx struct) - \param[in] hardness ice hardness - \param[in] gamma second invariant - \param[out] eta effective viscosity \f$ \eta \f$. - \param[out] derivative of the effective viscosity with respect - to the second invariant, \f$ \frac{\partial \eta}{\partial \gamma} \f$ - */ - PetscErrorCode (*viscosity)(void* ctx, PetscReal hardness, PetscReal gamma, - PetscReal *eta, PetscReal *deta); - - /*! Function evaluating the basal drag coefficient \f$ \tau_b \f$ as a - function of the basal yield stress \f$ \tau_c \f$ and - \f$ \gamma_b = u_b^2 + v_b^2 \f$. - - \param[in] ctx pointer to this context (BlatterQ1Ctx struct) - \param[in] tauc basal yield stress \f$ \tau_c \f$ - \param[in] gamma_b \f$ \gamma_b = u_b^2 + v_b^2 \f$ - \param[out] taub \f$ \tau_b \f$ - \param[out] dtaub derivative of \f$ \tau_b \f$ with respect to \f$ \gamma_b \f$ - */ - PetscErrorCode (*drag)(void *ctx, PetscReal tauc, PetscReal u, PetscReal v, - PetscReal *taub, PetscReal *dtaub); - } nonlinear; - - /*! 3D Q1 elements with the 2x2x2 (8-point) Gaussian quadrature */ - struct { - /*! Values of shape functions at quadrature points. - * - * chi[q][i] corresponds to quadrature point q, shape function i. - */ - PetscReal chi[8][8]; - /*! Derivatives of shape functions at quadrature points. - * - * dchi[q][i][d] corresponds to quadrature point q, shape function i, direction d. - */ - PetscReal dchi[8][8][3]; - } Q13D; - - /*! 2D Q1 elements with the 2x2 (4-point) Gaussian quadrature */ - struct { - /*! Values of shape functions at quadrature points. - * - * chi[q][i] corresponds to quadrature point q, shape function i. - */ - PetscReal chi[4][4]; - /*! Derivatives of shape functions at quadrature points. - * - * dchi[q][i][d] corresponds to quadrature point q, shape function i, direction d. - */ - PetscReal dchi[4][4][2]; - } Q12D; - - /*! Pointer to a PISM-side class that might contain parameters used - by viscosity() and drag(). */ - void *extra; -} BlatterQ1Ctx; - -PetscErrorCode BlatterQ1_begin_2D_parameter_access(DM da, PetscBool local, Vec *X_out, - PrmNode ***prm); -PetscErrorCode BlatterQ1_end_2D_parameter_access(DM da, PetscBool local, Vec *X_out, PrmNode ***prm); - -PetscErrorCode BlatterQ1_begin_hardness_access(DM da, PetscBool local, Vec *X_out, PetscScalar ****hardness); -PetscErrorCode BlatterQ1_end_hardness_access(DM da, PetscBool local, Vec *X_out, PetscScalar ****hardness); - -PetscErrorCode BlatterQ1_create(MPI_Comm com, DM pism_da, PetscInt Mz, - BlatterQ1Ctx *ctx, SNES *result); -#ifdef __cplusplus -} -#endif - - -#endif /* _BLATTER_IMPLEMENTATION_H_ */ diff --git a/src/stressbalance/blatter/CMakeLists.txt b/src/stressbalance/blatter/CMakeLists.txt index 58a6867726..fba78ba57a 100644 --- a/src/stressbalance/blatter/CMakeLists.txt +++ b/src/stressbalance/blatter/CMakeLists.txt @@ -1,6 +1,4 @@ add_library (blatter OBJECT - BlatterStressBalance.cc - Blatter_implementation.c grid_hierarchy.cc Poisson3.cc Blatter.cc diff --git a/src/stressbalance/blatter/FE3DTools.h b/src/stressbalance/blatter/FE3DTools.h deleted file mode 100644 index 0643fd22ab..0000000000 --- a/src/stressbalance/blatter/FE3DTools.h +++ /dev/null @@ -1,314 +0,0 @@ -/* Copyright (C) 2013, 2014 Jed Brown and the PISM Authors - * - * 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 3 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 _FE3DTools_H_ -#define _FE3DTools_H_ - -#include - -#if !defined __STDC_VERSION__ || __STDC_VERSION__ < 199901L -# if defined __cplusplus /* C++ restrict is nonstandard and compilers have inconsistent rules about where it can be used */ -# define restrict -# else -# define restrict PETSC_RESTRICT -# endif -#endif - -/*! \file FE3DTools.h - - This file contains helper macros and functions for setting up 3D - @f$ Q_1 @f$ finite element systems. - - In the finite element formulation, let @f$ \phi_i @f$ be a trial function, - i.e. an element of a basis for the space @f$ S_h @f$ approximating the - solution space and let @f$ \psi @f$ be a test function. - - Then we write - - \f[ - u = \sum_i U_i \phi_i. - \f] - - Trial functions @f$ \phi_i @f$ are piecewise-trilinear (linear in each of - @f$ x @f$ , @f$ y @f$ , @f$ z @f$ separately) on elements. - - Instead of defining @f$ \phi_i @f$ for each element, we define them once - on the **reference element**, the cube - @f$ [-1,1]\times[-1,1]\times[-1,1] @f$ and use an invertible map from this - element to a physical hexahedral element. - - **On** the reference element, pre-images of trial functions are - called @f$ \chi_i @f$ . These @f$ \chi_i @f$ are also known as *element basis - functions*, as opposed to *global* basis functions @f$ \phi_i @f$ . - - In the Galerkin formulation, test functions @f$ \psi @f$ are the same as - trial functions @f$ \phi_i @f$ , but it is still helpful to use different - letters for test and trial functions (this makes notation simpler). - */ - - -/*! Extract values of the 3D field x at nodes of an element (i,j,k) and store - them in the array n. - - Use this to *get* values of a field in a FEM assembly loop. -*/ -#define get_nodal_values_3d(x,i,j,k,n) do { \ - (n)[0] = (x)[i][j][k]; \ - (n)[1] = (x)[i+1][j][k]; \ - (n)[2] = (x)[i+1][j+1][k]; \ - (n)[3] = (x)[i][j+1][k]; \ - (n)[4] = (x)[i][j][k+1]; \ - (n)[5] = (x)[i+1][j][k+1]; \ - (n)[6] = (x)[i+1][j+1][k+1]; \ - (n)[7] = (x)[i][j+1][k+1]; \ - } while (0) - -/*! Extract memory addresses corresponding to nodes of an element (i,j,k) in - the field x and store them in the array n. - - Use this to *set* values of a field in a FEM assembly loop. - - This is used in the resudual evaluation code. -*/ -#define get_pointers_to_nodal_values_3d(x,i,j,k,n) do { \ - (n)[0] = &(x)[i][j][k]; \ - (n)[1] = &(x)[i+1][j][k]; \ - (n)[2] = &(x)[i+1][j+1][k]; \ - (n)[3] = &(x)[i][j+1][k]; \ - (n)[4] = &(x)[i][j][k+1]; \ - (n)[5] = &(x)[i+1][j][k+1]; \ - (n)[6] = &(x)[i+1][j+1][k+1]; \ - (n)[7] = &(x)[i][j+1][k+1]; \ - } while (0) - -/*! Extract values of a 2D field x at nodes of an element (i,j) and store them - in the array n. - - Use this to *get* values of a field in a FEM assembly loop. - - This is used in the code computing boundary integrals. -*/ -#define get_nodal_values_2d(x,i,j,n) do { \ - (n)[0] = (x)[i][j]; \ - (n)[1] = (x)[i+1][j]; \ - (n)[2] = (x)[i+1][j+1]; \ - (n)[3] = (x)[i][j+1]; \ - } while (0) - -/*! \brief Compute partial derivatives of @f$ z(\xi,\eta,\zeta) @f$ with - respect to @f$ \xi @f$ , @f$ \eta @f$ , and @f$ \zeta @f$ . - - The function @f$ z(\xi,\eta,\zeta) @f$ is defined by the map from the reference - element to a physical element, i.e. - \f[ - z (\xi, \eta, \zeta) = \sum_{j = 1}^8 z_j \cdot \chi_j (\xi,\eta,\zeta) . - \f] - - Here @f$ \xi @f$ , @f$ \eta @f$ and @f$ \zeta @f$ are range from @f$ -1 @f$ - to @f$ 1 @f$ , so that the triple @f$ (\xi,\eta,\zeta) @f$ describes an - arbitrary point in the reference element. - - These derivatives are used to compute the Jacobian of this map (they - appear in the third column). - - The output of this computation is used in compute_element_info() and nowhere - else. - - \param[in] dchi derivatives of element basis functions @f$ \phi @f$ with respect to - @f$ \xi @f$ , @f$ \eta @f$ , @f$ \zeta @f$ . - \param[in] zn[] z-coordinates of the nodes of the current element - \param[out] dz[] partial derivatives of z - */ -void compute_z_gradient(PetscReal dchi[][3], const PetscReal zn[], PetscReal dz[]) -{ - PetscInt i; - dz[0] = dz[1] = dz[2] = 0; - for (i = 0; i < 8; i++) { - dz[0] += dchi[i][0] * zn[i]; - dz[1] += dchi[i][1] * zn[i]; - dz[2] += dchi[i][2] * zn[i]; - } -} - -/*! Compute temporaries at a quadrature point. */ -/*! Compute the following: - - values of shape functions @f$ \phi_i @f$ at a given quadrature point - - values of partial derivatives of shape functions @f$ \frac{\partial \phi_i}{\partial x_j} @f$ - at a given quadrature point - - @f$ det(J)\cdot w @f$ factor (product of the determinant of the Jacobian of the map - from the reference element and the quadrature weight), at a given quadrature point. - (These are "modified" quadrature weights.) - - Let $J$ be the jacobian of the map from the reference element. - - This function computes @f$ J @f$ , @f$ det(J) @f$ , and @f$ J^{-1} @f$ . The determinant is then - multiplied by weights corresponding to the @f$ 2\times2\times2 @f$ Gaussian quadrature - (weights are hard-wired). - - The inverse @f$ J^{-1} @f$ is used to compute partial derivatives of - shape functions (with respect to @f$ x @f$ , @f$ y @f$ , @f$ z @f$ ) using - partial derivatives of reference element basis functions @f$ \chi_i @f$ (with - respect to @f$ \xi @f$ , @f$ \eta @f$ , @f$ \zeta @f$ ). - - In particular, @f$ \nabla \phi = J^{-1} (\nabla \chi) @f$ . - - **Note:** both the Jacobian and its inverse are stored *transposed* - here. (For no apparent reason.) - - \todo Quadrature weights are hard-wired! - - \param[in] chi values of element basis functions @f$ \chi_i @f$ at quadrature points - \param[in] dchi values of partial derivatives of @f$ \chi_i @f$ at quadrature points - \param[in] q index of the quadrature point - \param[in] dx,dy horizontal grid spacing (x- and y-direction) - \param[in] dz partial derivatives of z computed by calling compute_z_gradient() - \param[out] phi values of shape functions @f$ \phi_i @f$ at the quadrature point q - \param[out] dphi values of partial derivatives of shape functions @f$ \phi_i @f$ at q - \param[out] jw @f$ det(J)\cdot w @f$ - */ -void compute_element_info(PetscReal chi[8][8],PetscReal dchi[8][8][3], - PetscInt q, PetscReal dx, PetscReal dy, const PetscReal dz[restrict], - PetscReal phi[restrict], - PetscReal dphi[restrict][3], - PetscReal *restrict jw) -{ - const PetscReal jac[3][3] = {{dx / 2, 0, 0}, - {0, dy / 2, 0}, - {dz[0], dz[1], dz[2]}}, - ijac[3][3] = {{1 / jac[0][0], 0, 0}, - {0, 1 / jac[1][1], 0}, - { - jac[2][0] / (jac[0][0]*jac[2][2]), - jac[2][1] / (jac[1][1]*jac[2][2]), 1 / jac[2][2]}}, - jdet = jac[0][0]*jac[1][1]*jac[2][2]; - - for (int i = 0; i < 8; i++) { - const PetscReal *dchi_q = dchi[q][i]; - phi[i] = chi[q][i]; - dphi[i][0] = dchi_q[0]*ijac[0][0] + dchi_q[1]*ijac[1][0] + dchi_q[2]*ijac[2][0]; /* x-derivative */ - dphi[i][1] = dchi_q[0]*ijac[0][1] + dchi_q[1]*ijac[1][1] + dchi_q[2]*ijac[2][1]; /* y-derivative */ - dphi[i][2] = dchi_q[0]*ijac[0][2] + dchi_q[1]*ijac[1][2] + dchi_q[2]*ijac[2][2]; /* z-derivative */ - } - *jw = 1.0 * jdet; /* hard-wired quadrature weights */ -} - - -/*! Compute values of shape functions and their derivatives at quadrature points. - * - * Given coordinated of the nodes of the 2D @f$ Q_1 @f$ reference element - * \f{align*}{ - * \xi &= (-1, 1, 1, -1)\\ - * \eta &= (-1, -1, 1, 1) - * \f} - * we define 2D element basis functions - * \f{align*}{ - * \chi_i(\xi,\eta) &=\frac 1 4 (1 + \xi_i \xi)(1 + \eta_i \eta) - * \f} - * - * This function pre-computes values of these element basis functions - * and values of their derivatives at quadrature points for the - * @f$ 2\times2 @f$ Gaussian quadrature on the reference element. - * - * Note that with this choice of the reference element quadrature points are - * \f[ - * (\xi^*_i, \eta^*_i) = \left( \frac{\xi_i}{\sqrt3}, \frac{\eta_i}{\sqrt3} \right). - * \f] - * - * The partial derivatives are - * \f{align*}{ - * \frac{\partial \chi_i}{\partial \xi} &= \frac 1 4 \xi_i (1 + \eta_i \eta)\\ - * \frac{\partial \chi_i}{\partial \eta} &= \frac 1 4 \eta_i (1 + \xi_i \xi) - * \f} - * - * These 2D basis functions are used to approximate boundary integrals. - */ -void initialize_Q12D(PetscReal chi[4][4], PetscReal dchi[4][4][2]) -{ - /* Coordinates of the nodes of the reference element: */ - PetscReal xis[4] = {-1.0, 1.0, 1.0, -1.0}; - PetscReal etas[4] = {-1.0, -1.0, 1.0, 1.0}; - int q, j; - - for (q = 0; q < 4; ++q) { /* for all quadrature points... */ - /* compute coordinates of this quadrature point */ - PetscReal xi_q = xis[q]/sqrt(3.0), eta_q = etas[q]/sqrt(3.0); - - for (j = 0; j < 4; ++j) { - /* compute values of shape functions */ - chi[q][j] = 0.25 * (1.0 + xis[j]*xi_q) * (1.0 + etas[j]*eta_q); - - /* compute derivatives of shape functions */ - dchi[q][j][0] = 0.25 * xis[j] * (1.0 + etas[j] * eta_q); - dchi[q][j][1] = 0.25 * etas[j] * (1.0 + xis[j] * xi_q); - } - } -} - -/*! Compute values of shape functions and their derivatives at quadrature points. - * - * Given coordinated of the nodes of the 3D @f$ Q_1 @f$ reference element - * \f{align*}{ - * \xi &= (-1, 1, 1, -1, -1, 1, 1, -1)\\ - * \eta &= (-1, -1, 1, 1, -1, -1, 1, 1)\\ - * \zeta &= (-1, -1, -1, -1, 1, 1, 1, 1) - * \f} - * we define 3D element basis functions - * \f{align*}{ - * \chi_i(\xi,\eta,\zeta) &=\frac 1 8 (1 + \xi_i \xi)(1 + \eta_i \eta)(1 + \zeta_i \zeta) - * \f} - * - * This function pre-computes values of these element basis functions - * and values of their derivatives at quadrature points for the - * @f$ 2\times2\times2 @f$ Gaussian quadrature on the reference element. - * - * Note that with this choice of the reference element quadrature points are - * \f[ - * (\xi^*_i, \eta^*_i, \zeta^*_i) = \left( \frac{\xi_i}{\sqrt3}, \frac{\eta_i}{\sqrt3}, \frac{\zeta_i}{\sqrt3} \right). - * \f] - * - * The partial derivatives are - * \f{align*}{ - * \frac{\partial \chi_i}{\partial \xi} &= \frac 1 8 \xi_i (1 + \eta_i \eta)(1 + \zeta_i \zeta)\\ - * \frac{\partial \chi_i}{\partial \eta} &= \frac 1 8 \eta_i (1 + \xi_i \xi)(1 + \zeta_i \zeta)\\ - * \frac{\partial \chi_i}{\partial \zeta} &= \frac 1 8 \zeta_i (1 + \xi_i \xi)(1 + \eta_i \eta) - * \f} - */ -void initialize_Q13D(PetscReal chi[8][8], PetscReal dchi[8][8][3]) -{ - /* Coordinated of the nodes of the reference element: */ - PetscReal xis[8] = {-1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0}; - PetscReal etas[8] = {-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0}; - PetscReal zetas[8] = {-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0}; - int q, j; - - for (q = 0; q < 8; ++q) { /* for all quadrature points... */ - /* compute coordinates of this quadrature point */ - PetscReal xi_q = xis[q]/sqrt(3.0), eta_q = etas[q]/sqrt(3.0), zeta_q = zetas[q]/sqrt(3.0); - - for (j = 0; j < 8; ++j) { - /* compute values of shape functions */ - chi[q][j] = 0.125 * (1.0 + xis[j]*xi_q) * (1.0 + etas[j]*eta_q) * (1.0 + zetas[j]*zeta_q); - - /* compute derivatives of shape functions */ - dchi[q][j][0] = 0.125 * xis[j] * (1.0 + etas[j] * eta_q) * (1.0 + zetas[j] * zeta_q); - dchi[q][j][1] = 0.125 * etas[j] * (1.0 + xis[j] * xi_q) * (1.0 + zetas[j] * zeta_q); - dchi[q][j][2] = 0.125 * zetas[j] * (1.0 + xis[j] * xi_q) * (1.0 + etas[j] * eta_q); - } - } -} - -#endif /* _FE3DTools_H_ */