Port noise model

[chrisknewman]

We will use JL's method with the following:

1. tpetra_vector node_seed(x_overlapmap) -- node seeds at each timestep, communicate ghosts
2. need storage at elements for elem_seed. we have a concept of Teuchos::RCP elemmap; and Teuchos::RCP elemerror; in erroreestimator
3. need an array or std::vector for double gauss_val[NELEMENTS][NGAUSS]; to store random distributions. may need long long
4. need a view to push gauss_val to device. set gauss_val to all zero (or the view to all zero) if we don't need random numbers
5. need a bool from input file to determine if we need to compute randoms
6. need a method to populate gauss_val at each advance()
7. use bool to set up maps/vectors in constructor and to call the method
8. extract from view by elem_id, gp index, color. pass to resfuc. also pass volume to resfunc.
9. this should be in a class that has a method that provides gauss_val (it may look alot like error_estimator)

[chrisknewman]

include

include #define NELEMENTS 10// let's have 4 nodes/element, with consecutive pairs of elements sharing two nodes

define NNODES (NELEMENTS*2+2)

// 1----3----5--... // | | | // | | | // 0----2----4--...// just assume large number of Gauss points to enable verifications that they are uncorrelated

define NGAUSS 100000int main()

{ int initial_seed = 12345; // initialize node_seed array with random number between 1 and 10000 std::mt19937 gen(initial_seed); int node_seed[NNODES]; std::uniform_int_distribution<> udist(1,10000); for(int i=0;i<NNODES;i++) { node_seed[i]=udist(gen); //std::cout<<"val = "<<node_seed[i]<<std::endl; } // initialize one seed for each element int element_seed[NELEMENTS]; for(int i=0;i<NELEMENTS;i++) { element_seed[i]=0; // loop over 4 nodes of this element for(int node=0;node<4;node++) { int node_index=2i+node; element_seed[i]+=node_seed[node_index]; } std::cout<<"Element seed = "<<element_seed[i]<<std::endl; } // now set random values at Gauss points double gauss_val[NELEMENTS][NGAUSS]; for(int i=0;i<NELEMENTS;i++) { std::mt19937 mt(element_seed[i]); std::normal_distribution<> normal_dist(0,1); for(int ig=0;ig<NGAUSS;ig++) { gauss_val[i][ig]=normal_dist(mt); } } // verify suite of numbers in each element are uncorrelated double max_dot=0.; for(int i=0;i<NELEMENTS;i++) { for(int j=0;j<i;j++) { double dot=0.; for(int ig=0;ig<NGAUSS;ig++) { dot+=gauss_val[i][ig]gauss_val[j][ig]; } dot/=(double)NGAUSS; std::cout<<"Element "<<i<<","<<j<<": dot="<<dot<<std::endl; max_dot = dot>max_dot ? dot : max_dot; } } std::cout<<"Max. dot product: "<<max_dot<<std::endl; // compute self-correlation for comparison double ref_dot=0.; for(int ig=0;ig<NGAUSS;ig++) { ref_dot+=gauss_val[0][ig]*gauss_val[0][ig]; } ref_dot/=(double)NGAUSS; std::cout<<"Reference dot="<<ref_dot<<std::endl; return 0; }

[chrisknewman]

We could maybe skip the node_seed step and populate element_seed based on global id and timestep, then populate the random numbers at gauss points

[stvdwtt]

@chrisknewman, can we close this? Or is there cleaning-up you'd like to do?