SpatialOps
StencilBuilderCustom.cpp
1 /*
2  * The MIT License
3  *
4  * Copyright (c) 2016-2017 The University of Utah
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to
8  * deal in the Software without restriction, including without limitation the
9  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
10  * sell copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in
14  * all copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
21  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
22  * IN THE SOFTWARE.
23  */
24 
25 #include "StencilBuilder.h"
26 #include "FVStaggeredOperatorTypes.h"
27 #include "OneSidedOperatorTypes.h"
28 #include "FVStaggeredBCOp.h"
30 #include <spatialops/structured/Grid.h>
32 
33 namespace SpatialOps{
34 
35 #define REGISTER_BC_STENCILS( VOL ) \
36  opdb.register_new_bc_operator< BasicOpTypes<VOL>::InterpC2FX >(); \
37  opdb.register_new_bc_operator< BasicOpTypes<VOL>::InterpC2FY >(); \
38  opdb.register_new_bc_operator< BasicOpTypes<VOL>::InterpC2FZ >(); \
39  opdb.register_new_bc_operator< BasicOpTypes<VOL>::GradX >(); \
40  opdb.register_new_bc_operator< BasicOpTypes<VOL>::GradY >(); \
41  opdb.register_new_bc_operator< BasicOpTypes<VOL>::GradZ >(); \
42  opdb.register_new_bc_operator< BasicOpTypes<VOL>::DivX >(); \
43  opdb.register_new_bc_operator< BasicOpTypes<VOL>::DivY >(); \
44  opdb.register_new_bc_operator< BasicOpTypes<VOL>::DivZ >(); \
45  opdb.register_new_bc_operator< OperatorTypeBuilder<InterpolantX,VOL,VOL>::type >(); \
46  opdb.register_new_bc_operator< OperatorTypeBuilder<InterpolantY,VOL,VOL>::type >(); \
47  opdb.register_new_bc_operator< OperatorTypeBuilder<InterpolantZ,VOL,VOL>::type >(); \
48  opdb.register_new_bc_operator< OperatorTypeBuilder<GradientX, VOL,VOL>::type >(); \
49  opdb.register_new_bc_operator< OperatorTypeBuilder<GradientY, VOL,VOL>::type >(); \
50  opdb.register_new_bc_operator< OperatorTypeBuilder<GradientZ, VOL,VOL>::type >();
51 
52  //------------------------------------------------------------------
53 
54  template< typename FieldT, typename StencilT, typename CoefCollection >
55  void build_one_sided_stencils( OperatorDatabase& opdb,
56  const CoefCollection& coefs )
57  {
58  typedef typename OneSidedOpTypeBuilder<Gradient,StencilT,FieldT>::type GradOp;
59  opdb.register_new_operator( new GradOp( coefs ) );
60  }
61 
62  //------------------------------------------------------------------
63 
64  template< typename StencilT, typename CoefCollection >
65  void build_one_sided_stencils( OperatorDatabase& opdb,
66  const CoefCollection& coefs )
67  {
68  build_one_sided_stencils<SVolField,StencilT,CoefCollection>( opdb, coefs );
69  build_one_sided_stencils<XVolField,StencilT,CoefCollection>( opdb, coefs );
70  build_one_sided_stencils<YVolField,StencilT,CoefCollection>( opdb, coefs );
71  build_one_sided_stencils<ZVolField,StencilT,CoefCollection>( opdb, coefs );
72  }
73 
74  //------------------------------------------------------------------
75 
76  void build_custom_stencils( const unsigned int nx,
77  const unsigned int ny,
78  const unsigned int nz,
79  const double Lx,
80  const double Ly,
81  const double Lz,
82  OperatorDatabase& opdb )
83  {
84  const double dx = Lx/nx;
85  const double dy = Ly/ny;
86  const double dz = Lz/nz;
87 
88  //Coefficients:
89  NeboStencilCoefCollection<2> coefHalf = build_two_point_coef_collection( 0.5, 0.5 );
90  NeboStencilCoefCollection<2> coefDx = build_two_point_coef_collection( -1.0/dx, 1.0/dx );
91  NeboStencilCoefCollection<2> coefDy = build_two_point_coef_collection( -1.0/dy, 1.0/dy );
92  NeboStencilCoefCollection<2> coefDz = build_two_point_coef_collection( -1.0/dz, 1.0/dz );
93  NeboStencilCoefCollection<4> coefQuarter = build_four_point_coef_collection( 0.25, 0.25, 0.25, 0.25 );
94  NeboStencilCoefCollection<2> coefHalfDx = build_two_point_coef_collection( -0.5/dx, 0.5/dx );
95  NeboStencilCoefCollection<2> coefHalfDy = build_two_point_coef_collection( -0.5/dy, 0.5/dy );
96  NeboStencilCoefCollection<2> coefHalfDz = build_two_point_coef_collection( -0.5/dz, 0.5/dz );
97 
98  //___________________________________________________________________
99  // Finite Difference stencils:
100  //
101  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantX,SVolField,SVolField>::type( coefHalf ) );
102  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantY,SVolField,SVolField>::type( coefHalf ) );
103  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantZ,SVolField,SVolField>::type( coefHalf ) );
104  opdb.register_new_operator( new OperatorTypeBuilder<GradientX, SVolField,SVolField>::type( coefHalfDx ) );
105  opdb.register_new_operator( new OperatorTypeBuilder<GradientY, SVolField,SVolField>::type( coefHalfDy ) );
106  opdb.register_new_operator( new OperatorTypeBuilder<GradientZ, SVolField,SVolField>::type( coefHalfDz ) );
107 
108  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantX,XVolField,XVolField>::type( coefHalf ) );
109  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantY,XVolField,XVolField>::type( coefHalf ) );
110  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantZ,XVolField,XVolField>::type( coefHalf ) );
111  opdb.register_new_operator( new OperatorTypeBuilder<GradientX, XVolField,XVolField>::type( coefHalfDx ) );
112  opdb.register_new_operator( new OperatorTypeBuilder<GradientY, XVolField,XVolField>::type( coefHalfDy ) );
113  opdb.register_new_operator( new OperatorTypeBuilder<GradientZ, XVolField,XVolField>::type( coefHalfDz ) );
114 
115  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantX,YVolField,YVolField>::type( coefHalf ) );
116  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantY,YVolField,YVolField>::type( coefHalf ) );
117  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantZ,YVolField,YVolField>::type( coefHalf ) );
118  opdb.register_new_operator( new OperatorTypeBuilder<GradientX, YVolField,YVolField>::type( coefHalfDx ) );
119  opdb.register_new_operator( new OperatorTypeBuilder<GradientY, YVolField,YVolField>::type( coefHalfDy ) );
120  opdb.register_new_operator( new OperatorTypeBuilder<GradientZ, YVolField,YVolField>::type( coefHalfDz ) );
121 
122  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantX,ZVolField,ZVolField>::type( coefHalf ) );
123  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantY,ZVolField,ZVolField>::type( coefHalf ) );
124  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantZ,ZVolField,ZVolField>::type( coefHalf ) );
125  opdb.register_new_operator( new OperatorTypeBuilder<GradientX, ZVolField,ZVolField>::type( coefHalfDx ) );
126  opdb.register_new_operator( new OperatorTypeBuilder<GradientY, ZVolField,ZVolField>::type( coefHalfDy ) );
127  opdb.register_new_operator( new OperatorTypeBuilder<GradientZ, ZVolField,ZVolField>::type( coefHalfDz ) );
128 
129  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantX,VertexField,VertexField>::type( coefHalf ) );
130  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantY,VertexField,VertexField>::type( coefHalf ) );
131  opdb.register_new_operator( new OperatorTypeBuilder<InterpolantZ,VertexField,VertexField>::type( coefHalf ) );
132  opdb.register_new_operator( new OperatorTypeBuilder<GradientX, VertexField,VertexField>::type( coefHalfDx ) );
133  opdb.register_new_operator( new OperatorTypeBuilder<GradientY, VertexField,VertexField>::type( coefHalfDy ) );
134  opdb.register_new_operator( new OperatorTypeBuilder<GradientZ, VertexField,VertexField>::type( coefHalfDz ) );
135 
136  //___________________________________________________________________
137  // One sided stencils:
138  {
139  //One sided, 2 pt coefficients:
140  NeboStencilCoefCollection<2> coefHalfDx2Plus = build_two_point_coef_collection( -1.0/dx, 1.0/dx );
141  NeboStencilCoefCollection<2> coefHalfDx2Minus = build_two_point_coef_collection( 1.0/dx, -1.0/dx );
142  NeboStencilCoefCollection<2> coefHalfDy2Plus = build_two_point_coef_collection( -1.0/dy, 1.0/dy );
143  NeboStencilCoefCollection<2> coefHalfDy2Minus = build_two_point_coef_collection( 1.0/dy, -1.0/dy );
144  NeboStencilCoefCollection<2> coefHalfDz2Plus = build_two_point_coef_collection( -1.0/dz, 1.0/dz );
145  NeboStencilCoefCollection<2> coefHalfDz2Minus = build_two_point_coef_collection( 1.0/dz, -1.0/dz );
146 
147  //One sided, 3 pt coefficients:
148  NeboStencilCoefCollection<3> coefHalfDx3Plus = build_three_point_coef_collection( -1.5/dx, 2.0/dx, -0.5/dx );
149  NeboStencilCoefCollection<3> coefHalfDx3Minus = build_three_point_coef_collection( 1.5/dx, -2.0/dx, 0.5/dx );
150  NeboStencilCoefCollection<3> coefHalfDy3Plus = build_three_point_coef_collection( -1.5/dy, 2.0/dy, -0.5/dy );
151  NeboStencilCoefCollection<3> coefHalfDy3Minus = build_three_point_coef_collection( 1.5/dy, -2.0/dy, 0.5/dy );
152  NeboStencilCoefCollection<3> coefHalfDz3Plus = build_three_point_coef_collection( -1.5/dz, 2.0/dz, -0.5/dz );
153  NeboStencilCoefCollection<3> coefHalfDz3Minus = build_three_point_coef_collection( 1.5/dz, -2.0/dz, 0.5/dz );
154 
155  typedef UnitTriplet<XDIR>::type XPlus;
156  typedef UnitTriplet<YDIR>::type YPlus;
157  typedef UnitTriplet<ZDIR>::type ZPlus;
158  typedef XPlus::Negate XMinus;
159  typedef YPlus::Negate YMinus;
160  typedef ZPlus::Negate ZMinus;
161 
162  build_one_sided_stencils<OneSidedStencil2<XPlus > >( opdb, coefHalfDx2Plus );
163  build_one_sided_stencils<OneSidedStencil2<XMinus> >( opdb, coefHalfDx2Minus );
164  build_one_sided_stencils<OneSidedStencil2<YPlus > >( opdb, coefHalfDy2Plus );
165  build_one_sided_stencils<OneSidedStencil2<YMinus> >( opdb, coefHalfDy2Minus );
166  build_one_sided_stencils<OneSidedStencil2<ZPlus > >( opdb, coefHalfDz2Plus );
167  build_one_sided_stencils<OneSidedStencil2<ZMinus> >( opdb, coefHalfDz2Minus );
168 
169  build_one_sided_stencils<OneSidedStencil3<XPlus > >( opdb, coefHalfDx3Plus );
170  build_one_sided_stencils<OneSidedStencil3<XMinus> >( opdb, coefHalfDx3Minus );
171  build_one_sided_stencils<OneSidedStencil3<YPlus > >( opdb, coefHalfDy3Plus );
172  build_one_sided_stencils<OneSidedStencil3<YMinus> >( opdb, coefHalfDy3Minus );
173  build_one_sided_stencils<OneSidedStencil3<ZPlus > >( opdb, coefHalfDz3Plus );
174  build_one_sided_stencils<OneSidedStencil3<ZMinus> >( opdb, coefHalfDz3Minus );
175  }
176 
177  //---------------------------------------------------------------
178  // BC Stencils
179  REGISTER_BC_STENCILS( SVolField )
180  REGISTER_BC_STENCILS( XVolField )
181  REGISTER_BC_STENCILS( YVolField )
182  REGISTER_BC_STENCILS( ZVolField )
183  }
184 
185 } // namespace SpatialOps
SpatialField< YVol > YVolField
defines a volume field on the y-staggered volume
SpatialField< XVol > XVolField
defines a volume field on the x-staggered volume
SpatialField< SVol > SVolField
defines a volume field on the scalar volume.
SpatialField< ZVol > ZVolField
defines a volume field on the z-staggered volume