aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetInterp2/src/fasterp.hh
blob: d2a13a35fcfb4f360403bd13f52775d7a1e27a73 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
#ifndef FASTERP_HH
#define FASTERP_HH

#include <cstdlib>
#include <vector>

#include <cctk.h>

#include <defs.hh>
#include <dist.hh>
#include <typeprops.hh>
#include <vect.hh>



namespace CarpetInterp2 {
  
  using namespace std;
  
  
  
  int const dim = 3;
  
  // An interpolation point descriptor requires (3 * (max_order+1) +
  // 1) double precision values of memory
  int const max_order = 5;
  
  
  
  // Map C structures to MPI datatypes
  struct mpi_struct_descr_t {
    int          blocklength;
    MPI_Aint     displacement;
    MPI_Datatype type;
  };
  
  void
  create_mpi_datatype (size_t count,
                       mpi_struct_descr_t const descr[],
                       MPI_Datatype & newtype);
  
  
  
  // A global location, given by its global coordinates
  struct fasterp_glocs_t {
    vector<CCTK_REAL> coords[dim];
    fasterp_glocs_t (size_t const n)
    {
      for (int d=0; d<dim; ++d) {
        coords[d].resize(n);
      }
    }
    size_t size () const { return coords[0].size(); }
  };
  
  // A local location, given by map and local coordinates
  struct fasterp_llocs_t {
    vector<int> maps;
    vector<CCTK_REAL> coords[dim];
    fasterp_llocs_t (size_t const n)
    {
      maps.resize(n);
      for (int d=0; d<dim; ++d) {
        coords[d].resize(n);
      }
    }
    size_t size () const { return maps.size(); }
  };
  
  // An integer location, given by map, refinement level, and
  // component
  struct fasterp_iloc_t {
    int m, rl, c;
    
    // ivect idx;
    int ind3d;
    rvect offset;               // in terms of grid points
    
    static MPI_Datatype mpi_datatype ();
  };
  
  
  
  struct fasterp_dest_loc_t {
    // ivect idx;
    int ind3d;                  // destination grid point index
  };
  
  class fasterp_src_loc_t {
    CCTK_REAL coeffs[dim][max_order+1]; // interpolation coefficients
    bvect exact;
    
    // ivect idx;
    int ind3d;                  // source grid point offset
    
  public:
    void
    calc_stencil (fasterp_iloc_t const & iloc,
                  int order);
    void
    interpolate (ivect const & lsh,
                 int order,
                 vector<CCTK_REAL const *> const & varptrs,
                 CCTK_REAL * restrict vals)
      const;
    
  private:
    template <int O>
    void
    interpolate (ivect const & lsh,
                 vector<CCTK_REAL const *> const & varptrs,
                 CCTK_REAL * restrict vals)
      const;
    template <int O0, int O1, int O2>
    void
    interpolate (ivect const & lsh,
                 vector<CCTK_REAL const *> const & varptrs,
                 CCTK_REAL * restrict vals)
      const;
  };
  
  
  
  // A receive descriptor, describing what is received from other
  // processors
  struct recv_proc_t {
    int p;                      // sending processor
    int offset;
    int npoints;                // total number of received points
  };

  struct recv_descr_t {
    vector<recv_proc_t> procs;
    vector<int> procinds;
    int npoints;                // total number of received points
    
    vector<int> index;          // gather index list
  };
  
  // A send descriptor; describing what to send to other processors
  struct send_comp_t {
    // This structure does not exist for all components -- components
    // which are not accessed are not described, making this a sparse
    // data structure.  The field c contains the component number.
    vector<fasterp_src_loc_t> locs;
    int c;                      // source component
  };
  
  struct send_rl_t {
    vector<send_comp_t> comps;
    vector<int> compinds;
  };
  
  struct send_map_t {
    vector<send_rl_t> rls;
  };

  struct send_proc_t {
    // This structure does not exist for all processors -- processors
    // with which there is no communication are not described, making
    // this a sparse data structure.  The field p contains the
    // processor number.
    vector<send_map_t> maps;
    int p;                      // receiving processor
    int offset;
    int npoints;                // total number of sent points
  };
  
  struct send_descr_t {
    vector<send_proc_t> procs;
    // vector<int> procinds;
    int npoints;                // total number of sent points
  };
  
  
  
  class fasterp_setup_t {
    recv_descr_t recv_descr;
    send_descr_t send_descr;
    int order;
    
  public:
    fasterp_setup_t (cGH const * restrict cctkGH,
                     fasterp_glocs_t const & locations,
                     int order);
    
    ~ fasterp_setup_t ();
    
    void 
    interpolate (cGH const * restrict cctkGH,
                 vector<int> const & varinds,
                 vector<CCTK_REAL *> & values)
      const;
  };
  
  
  
} // namespace CarpetInterp2

#endif  // #define FASTERP_HH