Merge branch 'master' into 'master'

Refactor Science.cpp to split functions into separate files.

See merge request !1

src/Makefile

@@ -56,33 +56,42 @@ all: tests/test_deriv_x tests/test_write_x tests/test_esolve_x tests/test_heat_x
 
 .PHONY: clean
 clean:
-	rm -f *.o tests/*.o cases/*.o cases/*.src.? tests/*.src.?
+	rm -f *.o tests/*.o cases/*.o cases/*.src.? tests/*.src.? Science/*.o
+
+##
+## Short-names for source files
+##
 
 objfiles: $(shell ls *.cpp | sed -e 's/cpp/o/g')
 
+# SPINS files
+SPINS_BASE := TArray.o T_util.o Parformer.o ESolver.o Timestep.o NSIntegrator.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o Options.o BaseCase.o Sorter.o timing.o
+
+BASECASE_CPPS := $(wildcard BaseCase/*.cpp)
+BASECASE_OBJS := $(addprefix BaseCase/,$(notdir $(BASECASE_CPPS:.cpp=.o)))
+
+SCIENCE_CPPS := $(wildcard Science/*.cpp)
+SCIENCE_OBJS := $(addprefix Science/,$(notdir $(SCIENCE_CPPS:.cpp=.o)))
+
+##
+## Recipes
+##
+
 NSIntegrator.o: NSIntegrator.cpp NSIntegrator_impl.cc
 
 tests/test%.o: tests/test%.cpp
 	$(MPICXX) $(CFLAGS) -o $@ -c $<
 
-tests/test%_x: tests/test%.o tests/test%.src.o TArray.o Parformer.o T_util.o  ESolver.o Timestep.o NSIntegrator.o BaseCase.o Science.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o Options.o 
-	$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS) 
+tests/test%.x: tests/test%.o tests/test%.src.o ${SCIENCE_OBJS} ${SPINS_BASE}
+	$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
 
 cases/%.o: cases/%.cpp NSIntegrator_impl.cc NSIntegrator.hpp
 	$(MPICXX) $(CFLAGS) -o $@ -c  $<
 
-nonhydro_x: nonhydro_sw.o TArray.o T_util.o Parformer.o Splits.o Par_util.o Split_reader.o 
-	$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
-
-derek_x: derek.o TArray.o T_util.o Parformer.o ESolver.o Timestep.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o 
-	$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
-
-cases/%.x: cases/%.o cases/%.src.o TArray.o T_util.o Parformer.o ESolver.o Timestep.o NSIntegrator.o BaseCase.o Sorter.o \
-           Science.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o Options.o timing.o
+cases/%.x: cases/%.o cases/%.src.o ${SCIENCE_OBJS} ${SPINS_BASE}
 	$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
 
-cases/%_x: cases/%.o cases/%.src.o TArray.o T_util.o Parformer.o ESolver.o Timestep.o NSIntegrator.o BaseCase.o Sorter.o \
-           Science.o Splits.o Par_util.o Split_reader.o gmres.o gmres_1d_solver.o gmres_2d_solver.o grad.o multigrid.o Options.o timing.o
+cases/%_x: cases/%.o cases/%.src.o ${SCIENCE_OBJS} ${SPINS_BASE}
 	$(LD) $(LDFLAGS) -o $@ $^ $(LDLIBS)
 
 tests/test%.src.c: tests/test%.cpp CaseFileSource.c

src/Science.hpp

-/* WARNING: Science Content!
-
-   Collection of various analysis routines that are general enough to be useful
+   /* Collection of various analysis routines that are general enough to be useful
    over more than one project */
 
 #ifndef SCIENCE_HPP
@@ -11,6 +9,8 @@
 #include "NSIntegrator.hpp"
 
 
+// Global arrays to store quadrature weights
+extern Array<double,1> _quadw_x, _quadw_y, _quadw_z;
 
 // Marek's Overturning Diagnostic
 blitz::Array<double,3> overturning_2d(blitz::Array<double,3> const & rho, 
@@ -34,10 +34,12 @@ void compute_vort_z(TArrayn::DTArray & vortz, TArrayn::DTArray & u, TArrayn::DTA
 void compute_vorticity(TArrayn::DTArray & vortx, TArrayn::DTArray & vorty, TArrayn::DTArray & vortz,
         TArrayn::DTArray & u, TArrayn::DTArray & v, TArrayn::DTArray & w,
         TArrayn::Grad * gradient_op, const string * grid_type);
+
 // Enstrophy density
 void enstrophy_density(TArrayn::DTArray & enst, TArrayn::DTArray & u, TArrayn::DTArray & v,
         TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
         const int Nx, const int Ny, const int Nz);
+
 // Viscous dissipation
 void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray & v,
         TArrayn::DTArray & w, TArrayn::Grad * gradient_op, const string * grid_type,
@@ -47,6 +49,7 @@ void dissipation(TArrayn::DTArray & diss, TArrayn::DTArray & u, TArrayn::DTArray
 void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DTArray & quad3,
         int Nx, int Ny, int Nz, double Lx, double Ly, double Lz, double g, double rho_0, int iter,
         bool dimensional_rho = false, bool mapped = false, Array<double,1> hill = Array<double,1>());
+
 // Internal energy converted to BPE
 void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
         double kappa_rho, double rho_0, double g, int Nz, bool dimensional_rho = false);
@@ -56,12 +59,14 @@ void compute_quadweights(int szx, int szy, int szz,
       double Lx, double Ly, double Lz,
       NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
       NSIntegrator::DIMTYPE DIM_Z);
+
 const blitz::Array<double,1> * get_quad_x();
 const blitz::Array<double,1> * get_quad_y();
 const blitz::Array<double,1> * get_quad_z();
 
 // find which expansion to use based on field and boundary conditions
 void find_expansion(const string * grid_type, Transformer::S_EXP * expan, string deriv_filename);
+
 // switch trig function
 Transformer::S_EXP swap_trig( Transformer::S_EXP the_exp );
 
@@ -69,6 +74,7 @@ Transformer::S_EXP swap_trig( Transformer::S_EXP the_exp );
 void bottom_slope(TArrayn::DTArray & Hprime, TArrayn::DTArray & zgrid,
         TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
         const string * grid_type, const int Nx, const int Ny, const int Nz);
+
 // Top stresses
 void top_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & u,
         TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
@@ -76,6 +82,7 @@ void top_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & u,
 void top_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & v,
         TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
         const string * grid_type, const int Nz, const double visco);
+
 // Bottom stresses
 void bottom_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & Hprime,
         TArrayn::DTArray & u, TArrayn::DTArray & w, TArrayn::DTArray & temp,

src/Science/bottom_slope.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+
+// Bottom slope
+void bottom_slope(TArrayn::DTArray & Hprime, TArrayn::DTArray & zgrid,
+        TArrayn::DTArray & temp, TArrayn::Grad * gradient_op,
+        const string * grid_type, const int Nx, const int Ny, const int Nz) {
+    // Set-up
+    DTArray & z_x = *alloc_array(Nx,Ny,Nz);
+    blitz::Range all = blitz::Range::all();
+    blitz::firstIndex ii;
+    blitz::secondIndex jj;
+    blitz::thirdIndex kk;
+    S_EXP expan[3];
+    assert(gradient_op);
+
+    // get bottom topography
+    Array<double,1> xx(split_range(Nx));
+    xx = zgrid(all,0,0);
+    // put into temp array, and take derivative
+    temp = xx(ii) + 0*jj + 0*kk;
+    find_expansion(grid_type, expan, "zgrid");
+    gradient_op->setup_array(&temp,expan[0],expan[1],expan[2]);
+    gradient_op->get_dx(&z_x);
+    // flatten to get 2D array
+    Hprime(all,all,0) = z_x(all,all,0);
+    delete &z_x, xx;
+}

src/Science/bottom_stress_x.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+// Bottom Stress (along channel - x)
+void bottom_stress_x(TArrayn::DTArray & stress_x, TArrayn::DTArray & Hprime,
+        TArrayn::DTArray & u, TArrayn::DTArray & w, TArrayn::DTArray & temp,
+        TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
+        const double visco) {
+    // Set-up
+    blitz::Range all = blitz::Range::all();
+    S_EXP expan[3];
+    assert(gradient_op);
+    assert((grid_type[2] == "NO_SLIP") && "Surface stress requires NO_SLIP vertical (z) boundary condition");
+    // Along channel bottom stress can also be slightly inaccurate when x boundary condition is free slip
+    // since u or w do not necessarily satisfy Neumann BCs and the derivative has Gibbs at domain edges
+    // This is only true if there is velocity at left or right (in x) end wall corner, so is likely always small
+
+    if (mapped) {
+        // -du/dx
+        find_expansion(grid_type, expan, "u");
+        gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+        gradient_op->get_dx(&temp,false);
+        temp = (-1)*temp;
+        // dw/dz
+        find_expansion(grid_type, expan, "w");
+        gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+        gradient_op->get_dz(&temp,true);
+        // 2H'*(w_z-u_x)
+        stress_x(all,all,0) = 2*Hprime(all,all,0)*temp(all,all,0);
+
+        // dw/dx
+        find_expansion(grid_type, expan, "w");
+        gradient_op->setup_array(&w,expan[0],expan[1],expan[2]);
+        gradient_op->get_dx(&temp,false);
+        // du/dz
+        find_expansion(grid_type, expan, "u");
+        gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+        gradient_op->get_dz(&temp,true);
+        // (1-(H')^2)*(u_z+w_x)
+        stress_x(all,all,0) += (1-pow(Hprime(all,all,0),2))*temp(all,all,0);
+        // multiply by mu/(1+(H')^2)
+        stress_x = visco/(1+pow(Hprime,2))*stress_x;
+    } else {
+        // du/dz
+        find_expansion(grid_type, expan, "u");
+        gradient_op->setup_array(&u,expan[0],expan[1],expan[2]);
+        gradient_op->get_dz(&temp,false);
+        // bottom stress
+        stress_x(all,all,0) = visco*temp(all,all,0);
+    }
+}

src/Science/bottom_stress_y.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+// Bottom Stress (across channel - y)
+void bottom_stress_y(TArrayn::DTArray & stress_y, TArrayn::DTArray & Hprime,
+        TArrayn::DTArray & v, TArrayn::DTArray & temp,
+        TArrayn::Grad * gradient_op, const string * grid_type, const bool mapped,
+        const double visco) {
+    // Set-up
+    blitz::Range all = blitz::Range::all();
+    S_EXP expan[3];
+    assert(gradient_op);
+    assert((grid_type[2] == "NO_SLIP") && "Surface stress requires NO_SLIP vertical (z) boundary condition");
+    // Across channel bottom stress can also be slightly inaccurate when x boundary condition is free slip
+    // since v does not necessarily satisfy Neumann BCs and the derivative has Gibbs at domain edges
+    // This is only true if there is velocity at left or right (in x) end wall corner, so is likely always small
+
+    if (mapped) {
+        // dv/dx
+        find_expansion(grid_type, expan, "v");
+        gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+        gradient_op->get_dx(&temp,false);
+        // -v_x*H'
+        stress_y(all,all,0) = -temp(all,all,0)*Hprime(all,all,0);
+        // dv/dz
+        gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+        gradient_op->get_dz(&temp,false);
+        // add to -v_x*H'
+        stress_y(all,all,0) = temp(all,all,0) + stress_y(all,all,0);
+        // multiply by mu/sqrt(1+(H')^2)
+        stress_y = visco/pow(1+pow(Hprime,2),0.5)*stress_y;
+    } else {
+        // dv/dz
+        find_expansion(grid_type, expan, "v");
+        gradient_op->setup_array(&v,expan[0],expan[1],expan[2]);
+        gradient_op->get_dz(&temp,false);
+        // bottom stress
+        stress_y(all,all,0) = visco*temp(all,all,0);
+    }
+}

src/Science/compute_BPE_from_internal.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+// Internal energy converted to Background potential energy
+// See Winters et al. 1995 (Available potential energy and mixing in density-stratified fluids)
+// phi_i = - kappa * g * (integral of density at top boundary
+//                        minus the integral of density at bottom boundary)
+void compute_BPE_from_internal(double & phi_i, TArrayn::DTArray & rho,
+        double kappa_rho, double rho_0, double g, int Nz, bool dimensional_rho) {
+    // Tensor variables for indexing
+    blitz::firstIndex ii;
+    blitz::secondIndex jj;
+    blitz::Range all = blitz::Range::all();
+
+    if ( !dimensional_rho ) {
+        phi_i = -kappa_rho * g * pssum(sum(
+                    rho_0 * (rho(all,all,Nz-1) - rho(all,all,0))
+                    * (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
+    } else {
+        phi_i = -kappa_rho * g * pssum(sum(
+                    (rho(all,all,Nz-1) - rho(all,all,0))
+                    * (*get_quad_x())(ii) * (*get_quad_y())(jj) ));
+    }
+}

src/Science/compute_background_PE.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+bool compare_pairs( pair<double, double> a, pair<double, double> b ) {
+	return a.first < b.first;
+}
+
+// Compute Background Potential Energy (BPE)
+void compute_Background_PE(double & BPE_tot, TArrayn::DTArray & rho, TArrayn::DTArray & quad3,
+        int Nx, int Ny, int Nz, double Lx, double Ly, double Lz, double g,
+        double rho_0, int iter, bool dimensional_rho, bool mapped, Array<double,1> hill) {
+    // Tensor variables for indexing
+    blitz::firstIndex ii;
+    blitz::secondIndex jj;
+    blitz::thirdIndex kk;
+    // Set mpi parameters
+    int myrank, numprocs;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
+
+    // Arrays for sorting
+    static Array<double,3> sort_rho, sort_quad;
+    static vector<double> sort_hill(Nx), sort_dx(Nx);
+    static vector<double> Lx_partsum(Nx), hillvol_partsum(Nx);
+    double *local_vols =   new double[numprocs];
+    double *local_vols_r = new double[numprocs];
+    static double hill_max;
+    static double hill_vol;
+    static vector < pair<double, double> > height_width(Nx);
+
+    // Stuff to do once at the beginning
+    if (iter == 0) {
+        // adjust if mapped
+        if ( mapped ) {
+            // information about the hill
+            hill_vol = pssum(sum(hill*Ly*(*get_quad_x())));
+            hill_max = psmax(max(hill));
+
+            /* copy and sort the hill */
+            double *hill_tmp  = new double[Nx];
+            double *hill_tmp2 = new double[Nx];
+            double *dx_tmp  =   new double[Nx];
+            double *dx_tmp2 =   new double[Nx];
+            // create temporary empty arrays for holding the grid and hill
+            for (int II = 0; II < Nx; II++) {
+                if ( (II >= hill.lbound(firstDim)) and (II <= hill.ubound(firstDim)) ) {
+                    // populate these arrays with their processor's information
+                    hill_tmp2[II] = hill(II);
+                    dx_tmp2[II] = (*get_quad_x())(II);
+                } else {
+                    // else set to zero
+                    hill_tmp2[II] = 0.;
+                    dx_tmp2[II] = 0.;
+                }
+            }
+            // share across processors
+            MPI_Allreduce(hill_tmp2, hill_tmp, Nx, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+            MPI_Allreduce(dx_tmp2,   dx_tmp,   Nx, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+            // put into a double pair
+            for (int II = 0; II < Nx; II++) {
+                height_width[II].first  = hill_tmp[II];
+                height_width[II].second = dx_tmp[II];
+            }
+            // sort the height_width pair according to the height
+            sort(height_width.begin(), height_width.end(), compare_pairs);
+            // put sorted pairs into separate vectors
+            for (int II = 0; II < Nx; II++) {
+                sort_hill[II] = height_width[II].first;
+                sort_dx[II] = height_width[II].second;
+            }
+            // Compute cumulative sums of hill volume (reusing Lx_partsum)
+            for (int II = 0; II < Nx; II++) {
+                Lx_partsum[II] = sort_hill[II]*sort_dx[II]; 
+            }
+            std::partial_sum(Lx_partsum.begin(), Lx_partsum.end(), hillvol_partsum.begin());
+            // Compute cumulative sums of dxs
+            std::partial_sum(sort_dx.begin(), sort_dx.end(), Lx_partsum.begin());
+            // clean-up
+            delete[] dx_tmp;
+            delete[] dx_tmp2;
+            delete[] hill_tmp;
+            delete[] hill_tmp2;
+        }
+    }
+
+    // Compute sorted rho
+    sortarray(rho, quad3, sort_rho, sort_quad);
+
+    // Volume of memory-local portion of sorted array
+    double local_vol = sum(sort_quad);
+
+    // Share local volume information with other processors
+    // ie. processor 0 has lightest fluid
+    for (int II = 0; II < numprocs; II++) {
+        local_vols[II] = (II <= myrank) ? local_vol : 0;
+    }
+    MPI_Allreduce(local_vols, local_vols_r, numprocs, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+    // The total volume of fluid with greater or equal density
+    // to the lightest element on the local domain
+    double base_vol = local_vols_r[myrank];
+
+    // Find the depth at which base_vol will fill to
+    double tmpH, Area_star;
+    if ( !mapped ) {
+        tmpH = base_vol/Lx/Ly;
+    } else {
+        // check if base_vol will fill above the max hill height
+        if ((base_vol + hill_vol)/Lx/Ly >= hill_max) {
+            tmpH = (base_vol + hill_vol)/Lx/Ly;
+        } else {
+            // if it doesn't, find the depth where the fluid will fill to
+            int II = 0;
+            while ( (base_vol > Lx_partsum[II]*sort_hill[II] - hillvol_partsum[II]) and (II<Nx-1) ) {
+                II++;
+            }
+            assert(II>0 && "Something bad happened leading up to the tmpH calculation.");
+            // now subtract off the bit we went over by (draw yourself a picture, it works)
+            tmpH = sort_hill[II] - (sort_hill[II]*Lx_partsum[II] - base_vol - hillvol_partsum[II])/Lx_partsum[II-1];
+        }
+    }
+
+    // Now compute BPE from the sorted rho field
+    double dH = 0;
+    double BPE = 0;
+    int LL = Nx-1;
+    for (int II = sort_rho.lbound(firstDim); II <= sort_rho.ubound(firstDim); II++) {
+        for (int KK = sort_rho.lbound(thirdDim); KK <= sort_rho.ubound(thirdDim); KK++) {
+            for (int JJ = sort_rho.lbound(secondDim); JJ <= sort_rho.ubound(secondDim); JJ++) {
+                // Find the surface area of the water at depth tmpH
+                if ( mapped && (tmpH < hill_max) ) {
+                    while ( (sort_hill[LL] > tmpH) && (LL > 0) ) {
+                        LL--;
+                    }
+                    Area_star = Lx_partsum[LL]*Ly;
+                } else {
+                    Area_star = Lx*Ly;
+                }
+
+                // spread volume over domain and compute BPE for that cell
+                dH = sort_quad(II,JJ,KK)/Area_star;
+                if (dimensional_rho) {
+                    // rho is dimensional density
+                    BPE += g*sort_rho(II,JJ,KK)*(tmpH - 0.5*dH)*dH*Area_star;
+                } else {
+                    // rho is density anomaly
+                    BPE += g*rho_0*(1+sort_rho(II,JJ,KK))*(tmpH - 0.5*dH)*dH*Area_star;
+                }
+                tmpH -= dH;
+            }
+        }
+    }
+
+    // Share BPE with other processors
+    MPI_Allreduce(&BPE, &BPE_tot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+
+    // delete variables
+    delete[] local_vols;
+    delete[] local_vols_r;
+}

src/Science/compute_quadweights.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"
+#include "../Split_reader.hpp"
+#include "../T_util.hpp"
+#include "../Parformer.hpp"
+#include "../Sorter.hpp"
+#include <numeric>
+
+using blitz::Array;
+using blitz::cos;
+using namespace TArrayn;
+using namespace NSIntegrator;
+using namespace Transformer;
+
+// Global arrays to store quadrature weights
+Array<double,1> _quadw_x, _quadw_y, _quadw_z;
+
+// Compute quadrature weights
+void compute_quadweights(int szx, int szy, int szz, 
+      double Lx, double Ly, double Lz,
+      NSIntegrator::DIMTYPE DIM_X, NSIntegrator::DIMTYPE DIM_Y,
+      NSIntegrator::DIMTYPE DIM_Z) {
+   _quadw_x.resize(split_range(szx));
+   _quadw_y.resize(szy); _quadw_z.resize(szz);
+   if (DIM_X == NO_SLIP) {
+      blitz::firstIndex ii;
+      _quadw_x = 1;
+      for (int k = 1; k <= (szx-2)/2; k++) {
+         // From Trefethen, Spectral Methods in MATLAB
+         // clenshaw-curtis quadrature weights
+         _quadw_x -= 2*cos(2*k*M_PI*ii/(szx-1))/(4*k*k-1);
+      }
+      if ((szx%2))
+         _quadw_x -= cos(M_PI*ii)/((szx-1)*(szx-1)-1);
+      _quadw_x = 2*_quadw_x/(szx-1);
+      if (_quadw_x.lbound(firstDim) == 0) {
+         _quadw_x(0) = 1.0/(szx-1)/(szx-1);
+      }
+      if (_quadw_x.ubound(firstDim) == (szx-1)) {
+         _quadw_x(szx-1) = 1.0/(szx-1)/(szx-1);
+      }
+      _quadw_x *= Lx/2;
+   } else {
+      // Trapezoid rule
+      _quadw_x = Lx/szx;
+   }
+   if (DIM_Y == NO_SLIP) {
+      blitz::firstIndex ii;
+      _quadw_y = 1;
+      for (int k = 1; k <= (szy-2)/2; k++) {
+         // From Trefethen, Spectral Methods in MATLAB
+         // clenshaw-curtis quadrature weights
+         _quadw_y -= 2*cos(2*k*(M_PI*(ii)/(szy-1)))/(4*k*k-1);
+      }
+      if ((szy%2))
+         _quadw_y -= cos(M_PI*ii)/((szy-1)*(szy-1)-1);
+      _quadw_y = 2*_quadw_y/(szy-1);
+      _quadw_y(0) = 1.0/(szy-1)/(szy-1);
+      _quadw_y(szy-1) = 1.0/(szy-1)/(szy-1);
+      _quadw_y *= Ly/2;
+   } else {
+      // Trapezoid rule
+      _quadw_y = Ly/szy;
+   }
+   if (DIM_Z == NO_SLIP) {
+      blitz::firstIndex ii;
+      _quadw_z = 1;
+      for (int k = 1; k <= (szz-2)/2; k++) {
+         // From Trefethen, Spectral Methods in MATLAB
+         // clenshaw-curtis quadrature weights
+         _quadw_z -= 2*cos(2*k*(M_PI*(ii)/(szz-1)))/(4*k*k-1);
+      }
+      if ((szz%2))
+         _quadw_z -= cos(M_PI*ii)/((szz-1)*(szz-1)-1);
+      _quadw_z = 2*_quadw_z/(szz-1);
+      _quadw_z(0) = 1.0/(szz-1)/(szz-1);
+      _quadw_z(szz-1) = 1.0/(szz-1)/(szz-1);
+      _quadw_z *= Lz/2;
+   } else {
+      // Trapezoid rule
+      _quadw_z = Lz/szz;
+   }
+}

src/Science/compute_vort_x.cpp

+#include "../Science.hpp"
+#include "math.h"
+#include "../Par_util.hpp"
+#include "stdio.h"