tLingleClarkSerial.hh - 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
---
tLingleClarkSerial.hh (4438B)
---
1 // Copyright (C) 2007--2009, 2011, 2012, 2013, 2014, 2015, 2017, 2018, 2019 Ed Bueler 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 3 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 #ifndef LINGLECLARKSERIAL_H
20 #define LINGLECLARKSERIAL_H
21
22 #include <vector>
23
24 #include <petscvec.h>
25 #include <fftw3.h>
26 #include <vector>
27
28 #include "pism/util/petscwrappers/Vec.hh"
29 #include "pism/util/Logger.hh"
30
31 namespace pism {
32
33 class Config;
34
35 namespace bed {
36
37 //! Class implementing the bed deformation model described in [@ref BLKfastearth].
38 /*!
39 This class implements the [@ref LingleClark] bed deformation model by a Fourier
40 spectral collocation method, as described in [@ref BLKfastearth]. (The former
41 reference is where the continuum model arose, and a flow-line application is given.
42 The latter reference describes a new, fast method and gives verification results.
43 See also [@ref BLK2006earth] if more technical detail and/or Matlab programs are desired.)
44
45 Both a viscous half-space model (with elastic
46 lithosphere) and a spherical elastic model are computed. They are superposed
47 because the underlying earth model is linear.
48
49 The class assumes that the supplied Petsc Vecs are *sequential*. It is expected to be
50 run only on processor zero (or possibly by each processor once each processor
51 owns the entire 2D gridded load thicknesses and bed elevations.)
52
53 This model always assumes that we start with no load. Note that this does not mean that we
54 starting state is the equilibrium: the viscous plate may be "pre-bent" by using a provided
55 displacement field or by computing its displacement using an uplift field.
56 */
57 class LingleClarkSerial {
58 public:
59 LingleClarkSerial(Logger::ConstPtr log,
60 const Config &config,
61 bool include_elastic,
62 int Mx, int My,
63 double dx, double dy,
64 int Nx, int Ny);
65 ~LingleClarkSerial();
66
67 void init(Vec viscous_displacement,
68 Vec elastic_displacement);
69
70 void bootstrap(Vec thickness, Vec uplift);
71
72 void step(double dt_seconds, Vec H);
73
74 Vec total_displacement() const;
75
76 Vec viscous_displacement() const;
77
78 Vec elastic_displacement() const;
79
80 void compute_load_response_matrix(fftw_complex *output);
81 private:
82 void compute_elastic_response(Vec H, Vec dE);
83
84 void uplift_problem(Vec load_thickness, Vec bed_uplift, Vec output);
85
86 void precompute_coefficients();
87
88 void update_displacement(Vec V, Vec dE, Vec dU);
89
90 bool m_include_elastic;
91 // grid size
92 int m_Mx;
93 int m_My;
94 // grid spacing
95 double m_dx;
96 double m_dy;
97 //! load density (for computing load from its thickness)
98 double m_load_density;
99 //! mantle density
100 double m_mantle_density;
101 //! mantle viscosity
102 double m_eta;
103 //! lithosphere flexural rigidity
104 double m_D;
105
106 // acceleration due to gravity
107 double m_standard_gravity;
108
109 // size of the extended grid
110 int m_Nx;
111 int m_Ny;
112
113 // indices into extended grid for the corner of the physical grid
114 int m_i0_offset;
115 int m_j0_offset;
116
117 // half-lengths of the extended (FFT, spectral) computational domain
118 double m_Lx;
119 double m_Ly;
120
121 // Coefficients of derivatives in Fourier space
122 std::vector<double> m_cx, m_cy;
123
124 // viscous displacement on the extended grid
125 petsc::Vec m_Uv;
126
127 // elastic plate displacement
128 petsc::Vec m_Ue;
129
130 // total (viscous and elastic) plate displacement
131 petsc::Vec m_U;
132
133 fftw_complex *m_fftw_input;
134 fftw_complex *m_fftw_output;
135 fftw_complex *m_loadhat;
136 fftw_complex *m_lrm_hat;
137
138 fftw_plan m_dft_forward;
139 fftw_plan m_dft_inverse;
140
141 void tweak(Vec load_thickness, Vec U, int Nx, int Ny, double time);
142
143 Logger::ConstPtr m_log;
144 };
145
146 } // end of namespace bed
147 } // end of namespace pism
148
149 #endif /* LINGLECLARKSERIAL_H */