SpatialOps
heat_equation.cpp
Go to the documentation of this file.
1 
54 #include <spatialops/structured/Grid.h>
56 
57 using namespace SpatialOps;
58 
59 // If we are compiling with GPU CUDA support, create fields on the device.
60 // Otherwise, create them on the host.
61 #ifdef ENABLE_CUDA
62 # define LOCATION GPU_INDEX
63 #else
64 # define LOCATION CPU_INDEX
65 #endif
66 
67 //==============================================================================
68 
69 template<typename FieldT>
70 void
71 initialize_thermal_diffusivity( const Grid& grid, FieldT & alpha )
72 {
73  SpatFldPtr<FieldT> x = SpatialFieldStore::get<FieldT>( alpha );
74  SpatFldPtr<FieldT> y = SpatialFieldStore::get<FieldT>( alpha );
75 
76  grid.set_coord<XDIR>( *x );
77  grid.set_coord<YDIR>( *y );
78  //alpha <<= 0.1 + 0.4 * (*x + *y + 1.0);
79  //alpha <<= 0.1 + 0.4 * (*x + *y);
80  //alpha <<= (*x + *y + .2) /2;
81  alpha <<= 1.0;
82 
83 # ifdef ENABLE_CUDA
84  alpha.add_device( CPU_INDEX ); // transfer to facilitate printing its values
85 # endif
86  std::cout << "Initial alpha:" << std::endl;
87  //
88  //print_field( alpha, std::cout, true );
89 }
90 
91 //==============================================================================
92 
93 void
94 initialize_mask_points( const Grid& grid,
95  std::vector<IntVec>& leftSet,
96  std::vector<IntVec>& rightSet,
97  std::vector<IntVec>& topSet,
98  std::vector<IntVec>& bottomSet )
99 {
100  for( int i=-1; i<=grid.extent(1); ++i ){
101  leftSet .push_back( IntVec(-1, i, 0) );
102  rightSet.push_back( IntVec(grid.extent(0), i, 0) );
103  }
104 
105  for(int i = -1; i <= grid.extent(0); i++){
106  topSet .push_back( IntVec(i, -1, 0) );
107  bottomSet.push_back( IntVec(i, grid.extent(1), 0) );
108  }
109 }
110 
111 //==============================================================================
112 
113 template<typename FieldT>
114 double
115 find_deltaT( const FieldT& alpha, const Grid& grid )
116 {
117  const double deltaX = grid.spacing<XDIR>();
118  const double deltaY = grid.spacing<YDIR>();
119  const double sqrdDeltaX = deltaX * deltaX;
120  const double sqrdDeltaY = deltaY * deltaY;
121  const double sqrdDeltaXYmult = sqrdDeltaX * sqrdDeltaY;
122  const double sqrdDeltaXYplus = sqrdDeltaX + sqrdDeltaY;
123 
124  return 0.25 * sqrdDeltaXYmult / ( sqrdDeltaXYplus * nebo_max(alpha) );
125 }
126 
127 //==============================================================================
128 
129 int main()
130 {
131  typedef SVolField FieldT;
132 
133  //----------------------------------------------------------------------------
134  // Define the domain size and number of points
135  const DoubleVec domainLength(1,1,1); // a cube of unit length
136  const IntVec fieldDim( 6, 6, 1 ); // 6 x 6 x 1 points
137 
138  //----------------------------------------------------------------------------
139  // Create fields
140  const GhostData nghost(1);
141  const BoundaryCellInfo bcInfo = BoundaryCellInfo::build<FieldT>( true, true, true );
142  const MemoryWindow window( get_window_with_ghost( fieldDim, nghost, bcInfo) );
143 
144  FieldT phi( window, bcInfo, nghost, NULL, InternalStorage, LOCATION );
145  FieldT rhs( window, bcInfo, nghost, NULL, InternalStorage, LOCATION );
146  FieldT alpha( window, bcInfo, nghost, NULL, InternalStorage, LOCATION );
147 
148  const Grid grid( fieldDim, domainLength );
149 
150  //----------------------------------------------------------------------------
151  // Initialize alpha, thermal diffusivity
152  initialize_thermal_diffusivity( grid, alpha );
153 
154  //----------------------------------------------------------------------------
155  // Build and initialize masks:
156  std::vector<IntVec> leftSet, rightSet, topSet, bottomSet;
157 
158  initialize_mask_points( grid, leftSet, rightSet, topSet, bottomSet );
159 
160  SpatialMask<FieldT> left ( phi, leftSet );
161  SpatialMask<FieldT> right ( phi, rightSet );
162  SpatialMask<FieldT> top ( phi, topSet );
163  SpatialMask<FieldT> bottom( phi, bottomSet );
164 
165 # ifdef ENABLE_CUDA
166  // Masks are created on CPU so we need to explicitly transfer them to GPU
167  left .add_consumer( GPU_INDEX );
168  right .add_consumer( GPU_INDEX );
169  top .add_consumer( GPU_INDEX );
170  bottom.add_consumer( GPU_INDEX );
171 # endif
172 
173  //----------------------------------------------------------------------------
174  // Build stencils:
175  OperatorDatabase opDB; // holds stencils that can be retrieved easily
176  build_stencils( grid, opDB ); // builds stencils and places them in opDB
177 
178  typedef BasicOpTypes<SVolField>::DivX DivX; // x-divergence operator type
179  typedef BasicOpTypes<SVolField>::GradX GradX; // x-gradient operator type
180  typedef BasicOpTypes<SVolField>::InterpC2FX InterpX; // x-interpolant operator type
181 
182  const DivX& divX = *opDB.retrieve_operator<DivX >(); // retrieve the DivX operator
183  const GradX& gradX = *opDB.retrieve_operator<GradX >(); // retrieve the GradX operator
184  const InterpX& interpX = *opDB.retrieve_operator<InterpX>(); // retrieve the InterpX operator
185 
186  typedef BasicOpTypes<SVolField>::DivY DivY; // y-divergence operator type
187  typedef BasicOpTypes<SVolField>::GradY GradY; // y-gradient operator type
188  typedef BasicOpTypes<SVolField>::InterpC2FY InterpY; // y-interpolant operator type
189 
190  const DivY& divY = *opDB.retrieve_operator<DivY >(); // retrieve the DivY operator
191  const GradY& gradY = *opDB.retrieve_operator<GradY >(); // retrieve the GradY operator
192  const InterpY& interpY = *opDB.retrieve_operator<InterpY>(); // retrieve the InterpY operator
193 
194  //----------------------------------------------------------------------------
195  // Determine a safe deltaT:
196  const double deltaT = find_deltaT( alpha, grid );
197 
198  //----------------------------------------------------------------------------
199  // Initialize phi:
200  phi <<= cond( left, 10.0 )
201  ( right, 0.0 )
202  ( 5.0 );
203 
204  std::cout << "Initial phi:" << std::endl;
205  print_field( phi, std::cout, true );
206 
207  //----------------------------------------------------------------------------
208  // Take time steps:
209  const int nSteps = 40;
210  for( int i=0; i<=nSteps; ++i ){
211 
212  // update the solution using a forward-Euler update
213  phi <<= phi + deltaT * (
214  divX( interpX(alpha) * gradX(phi) )
215  + divY( interpY(alpha) * gradY(phi) ) );
216 
217  // Reset boundaries:
218  phi <<= cond( left, 10.0 )
219  ( right, 0.0 )
220  ( top || bottom, 5.0 )
221  ( phi );
222 
223  // print current state:
224  // print_field( phi, std::cout, true );
225  }
226 
227  return 0;
228 }
double spacing() const
obtain the grid spacing in the requested direction
Definition: Grid.cpp:53
Provides tools to index into a sub-block of memory.
Definition: MemoryWindow.h:63
MemoryWindow get_window_with_ghost(const IntVec &localDim, const GhostData &ghost, const BoundaryCellInfo &bc)
Obtain the memory window for a field on a patch that is a single, contiguous memory block...
Definition: MemoryWindow.h:314
void print_field(const Field &f, std::ostream &os, const bool addFormat=false)
print the values of a field (and ghost cells) to standard output
Definition: FieldHelper.h:281
void set_coord(FieldT &f) const
set the coordinates on the given field.
Definition: Grid.cpp:61
void build_stencils(const unsigned int nx, const unsigned int ny, const unsigned int nz, const double Lx, const double Ly, const double Lz, OperatorDatabase &opdb)
builds commonly used stencil operators
Abstracts a mask.
Definition: SpatialMask.h:70
OpT * retrieve_operator(const int id=-1) const
Retrieve an operator of the requested type.
Provides a simple interface to set coordinate fields.
Definition: Grid.h:36
Provides typedefs for common operator types on a given volume.
Holds information about the number of ghost cells on each side of the domain.
Definition: GhostData.h:54
Defines a type for the x-direction.
Defines a type for the y-direction.
Provides a database to hold operators of any type.
Abstracts a field.
Definition: SpatialField.h:132
Provides information about boundary cells for various fields.
Wrapper for pointers to SpatialField objects. Provides reference counting.