tfracture_density.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tfracture_density.cc (2829B)
---
1 // Copyright (C) 2011-2019 Torsten Albrecht and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 2 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 #include "IceModel.hh"
20
21 #include "pism/energy/EnergyModel.hh"
22 #include "pism/stressbalance/ShallowStressBalance.hh"
23 #include "pism/fracturedensity/FractureDensity.hh"
24
25 namespace pism {
26
27 void IceModel::update_fracture_density() {
28 // generate the BC mask for the fracture density model
29 //
30 // This mask contains ones at the in-flow boundary according to the SSA Dirichlet BC
31 // mask (if it is used) and at grounded grid points if
32 // fracture_density.include_grounded_ice is not set.
33 auto &bc_mask = m_work2d[0];
34 {
35 bool do_fracground = m_config->get_flag("fracture_density.include_grounded_ice");
36 const bool dirichlet_bc = m_config->get_flag("stress_balance.ssa.dirichlet_bc");
37
38 bc_mask.set(0.0);
39
40 IceModelVec::AccessList list{&bc_mask, &m_geometry.cell_type};
41
42 if (dirichlet_bc) {
43 list.add(m_ssa_dirichlet_bc_mask);
44 list.add(m_ssa_dirichlet_bc_values);
45 }
46
47 for (Points p(*m_grid); p; p.next()) {
48 const int i = p.i(), j = p.j();
49
50 if (m_geometry.cell_type.grounded(i, j) and not do_fracground) {
51 bc_mask(i, j) = 1.0;
52 }
53
54 if (dirichlet_bc) {
55 if (m_ssa_dirichlet_bc_mask(i, j) > 0.5 and
56 (m_ssa_dirichlet_bc_values(i, j).u != 0.0 or
57 m_ssa_dirichlet_bc_values(i, j).v != 0.0)) {
58 bc_mask(i, j) = 1.0;
59 }
60 }
61 } // end of the loop over grid points
62 }
63
64 // compute the vertically-averaged ice hardness
65 auto &hardness = m_work2d[1];
66 {
67 rheology::averaged_hardness_vec(*m_stress_balance->shallow()->flow_law(),
68 m_geometry.ice_thickness,
69 m_energy_model->enthalpy(),
70 hardness);
71 }
72
73 // This model has the same time-step restriction as the mass transport code so we don't
74 // check if this time step is short enough.
75 m_fracture->update(m_dt, m_geometry,
76 m_stress_balance->shallow()->velocity(),
77 hardness, bc_mask);
78 }
79
80 } // end of namespace pism