Merge branch 'master' of github.com:JamesEMcClure/LBPM-WIA
This commit is contained in:
JamesEMcclure
commit
5b6628a5c1
|
|
@ -90,7 +90,7 @@ CHECK_ENABLE_FLAG( USE_DOXYGEN 1 )
|
|||
CHECK_ENABLE_FLAG( USE_LATEX 1 )
|
||||
FILE( MAKE_DIRECTORY "${${PROJ}_INSTALL_DIR}/doc" )
|
||||
IF ( USE_DOXYGEN )
|
||||
SET( DOXYFILE_LATEX YES )
|
||||
SET( DOXYFILE_LATEX NO )
|
||||
SET( DOXYFILE_IN "${${PROJ}_SOURCE_DIR}/doxygen/Doxyfile.in" )
|
||||
SET( DOXY_HEADER_FILE "${${PROJ}_SOURCE_DIR}/doxygen/html/header.html" )
|
||||
SET( DOXY_FOOTER_FILE "${${PROJ}_SOURCE_DIR}/doxygen/html/footer.html" )
|
||||
|
|
|
|||
|
|
@ -102,7 +102,7 @@ std::shared_ptr<IO::Variable> getVariable( const std::string &path, const std::s
|
|||
* @brief Reformat the variable to match the mesh
|
||||
* @details This function modifies the dimensions of the array to match the mesh
|
||||
* @param[in] mesh The underlying mesh
|
||||
* @param[in/out] variable The variable name to read
|
||||
* @param[in,out] var The variable name to read
|
||||
*/
|
||||
void reformatVariable( const IO::Mesh &mesh, IO::Variable &var );
|
||||
|
||||
|
|
|
|||
|
|
@ -105,6 +105,7 @@ public:
|
|||
* @param[in] type The element type
|
||||
* @param[in] NumElements The number of elements
|
||||
* @param[in] dofMap The connectivity information (type x NumElements)
|
||||
* @param[in] NumNodes The number of nodes
|
||||
* @param[in] x The x coordinates or the xy/xyz coordinates
|
||||
* @param[in] y The y coordinates (may be null)
|
||||
* @param[in] z The z coordinates (may be null)
|
||||
|
|
|
|||
|
|
@ -138,6 +138,8 @@ Array<TYPE> getAtt( int fid, const std::string &att );
|
|||
* @brief Write the dimensions
|
||||
* @details This function writes the grid dimensions to netcdf.
|
||||
* @param fid Handle to the open file
|
||||
* @param names
|
||||
* @param dims
|
||||
*/
|
||||
std::vector<int> defDim(
|
||||
int fid, const std::vector<std::string> &names, const std::vector<int> &dims );
|
||||
|
|
@ -147,6 +149,10 @@ std::vector<int> defDim(
|
|||
* @brief Write a variable
|
||||
* @details This function writes a variable to netcdf.
|
||||
* @param fid Handle to the open file
|
||||
* @param var Variable to read
|
||||
* @param dimids
|
||||
* @param data
|
||||
* @param rank_info
|
||||
*/
|
||||
template<class TYPE>
|
||||
void write( int fid, const std::string &var, const std::vector<int> &dimids,
|
||||
|
|
|
|||
|
|
@ -64,7 +64,7 @@ void close( DBfile *fid );
|
|||
* @param[in] fid Handle to the open file
|
||||
* @param[in] name Name of variable
|
||||
*/
|
||||
DataType varDataType( DBfile *dbfile, const std::string &name );
|
||||
DataType varDataType( DBfile *fid, const std::string &name );
|
||||
|
||||
|
||||
/*!
|
||||
|
|
@ -265,8 +265,6 @@ void writeMultiMesh( DBfile *fid, const std::string &meshname,
|
|||
* @param[in] varname Mesh name
|
||||
* @param[in] subVarNames Names of the sub meshes in the form "filename:meshname"
|
||||
* @param[in] subVarTypes Type of each submesh
|
||||
* @param[in] ndim Dimension of variable (used to determine suffix)
|
||||
* @param[in] nvar Number of subvariables (used to determine suffix)
|
||||
*/
|
||||
void writeMultiVar( DBfile *fid, const std::string &varname,
|
||||
const std::vector<std::string> &subVarNames, const std::vector<int> &subVarTypes );
|
||||
|
|
|
|||
|
|
@ -207,6 +207,7 @@ double FlowAdaptor::UpdateFractionalFlow(ScaLBL_ColorModel &M){
|
|||
|
||||
void FlowAdaptor::Flatten(ScaLBL_ColorModel &M){
|
||||
|
||||
ScaLBL_D3Q19_Init(M.fq, M.Np);
|
||||
ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, 0, M.ScaLBL_Comm->LastExterior(), M.Np);
|
||||
ScaLBL_PhaseField_Init(M.dvcMap, M.Phi, M.Den, M.Aq, M.Bq, M.ScaLBL_Comm->FirstInterior(), M.ScaLBL_Comm->LastInterior(), M.Np);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -38,8 +38,8 @@ typedef Array<BlobIDType> BlobIDArray;
|
|||
* @param[in] SignDist SignDist
|
||||
* @param[in] vF vF
|
||||
* @param[in] vS vS
|
||||
* @param[in] S S
|
||||
* @param[out] LocalBlobID The ids of the blobs
|
||||
* @param[in] periodic Optional value
|
||||
* @return Returns the number of blobs
|
||||
*/
|
||||
int ComputeLocalBlobIDs( const DoubleArray& Phase, const DoubleArray& SignDist,
|
||||
|
|
@ -64,12 +64,13 @@ int ComputeLocalPhaseComponent( const IntArray &PhaseID, int &VALUE, IntArray &C
|
|||
* @param[in] nx Number of elements in the x-direction
|
||||
* @param[in] ny Number of elements in the y-direction
|
||||
* @param[in] nz Number of elements in the z-direction
|
||||
* @param[in] rank_info MPI communication info
|
||||
* @param[in] Phase Phase
|
||||
* @param[in] SignDist SignDist
|
||||
* @param[in] vF vF
|
||||
* @param[in] vS vS
|
||||
* @param[in] S S
|
||||
* @param[out] LocalBlobID The ids of the blobs
|
||||
* @param[out] GlobalBlobID The ids of the blobs
|
||||
* @param[in] comm MPI communicator
|
||||
* @return Returns the number of blobs
|
||||
*/
|
||||
int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_info,
|
||||
|
|
@ -84,10 +85,11 @@ int ComputeGlobalBlobIDs( int nx, int ny, int nz, const RankInfoStruct& rank_inf
|
|||
* @param[in] nx Number of elements in the x-direction
|
||||
* @param[in] ny Number of elements in the y-direction
|
||||
* @param[in] nz Number of elements in the z-direction
|
||||
* @param[in] rank_in MPI communication info
|
||||
* @param[in] rank_info MPI communication info
|
||||
* @param[in] PhaseID Array that identifies the phases
|
||||
* @param[in] VALUE Identifier for the phase to decompose
|
||||
* @param[out] GlobalBlobID The ids of the blobs for the phase
|
||||
* @param[in] comm The communicator to use
|
||||
* @return Return the number of components in the specified phase
|
||||
*/
|
||||
int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& rank_info,
|
||||
|
|
@ -98,10 +100,8 @@ int ComputeGlobalPhaseComponent( int nx, int ny, int nz, const RankInfoStruct& r
|
|||
* @brief Reorder the blobs
|
||||
* @details Reorder the blobs based on the number of cells they contain
|
||||
* largest first.
|
||||
* @param[in] nx Number of elements in the x-direction
|
||||
* @param[in] ny Number of elements in the y-direction
|
||||
* @param[in] nz Number of elements in the z-direction
|
||||
* @param[in/out] ID The ids of the blobs
|
||||
* @param[in,out] ID The ids of the blobs
|
||||
* @param[in] comm MPI communicator
|
||||
*/
|
||||
void ReorderBlobIDs( BlobIDArray& ID, const Utilities::MPI& comm );
|
||||
|
||||
|
|
@ -133,8 +133,12 @@ struct ID_map_struct {
|
|||
* @details This functions computes the map of blob ids between iterations
|
||||
* @return Returns the map of the blob ids. Each final blob may have no source
|
||||
* ids, one parent, or multiple parents. Each src id may be a parent for multiple blobs.
|
||||
* @param[in] nx Number of elements in the x-direction
|
||||
* @param[in] ny Number of elements in the y-direction
|
||||
* @param[in] nz Number of elements in the z-direction
|
||||
* @param[in] ID1 The blob ids at the first timestep
|
||||
* @param[in] ID2 The blob ids at the second timestep
|
||||
* @param[in] comm The communicator to use
|
||||
*/
|
||||
ID_map_struct computeIDMap( int nx, int ny, int nz, const BlobIDArray& ID1, const BlobIDArray& ID2, const Utilities::MPI& comm );
|
||||
|
||||
|
|
@ -143,7 +147,7 @@ ID_map_struct computeIDMap( int nx, int ny, int nz, const BlobIDArray& ID1, cons
|
|||
* @brief Compute the new global ids based on the map
|
||||
* @details This functions computes the time-consistent global ids for the
|
||||
* current global id index
|
||||
* @param[in/out] map The timestep mapping for the ids
|
||||
* @param[in,out] map The timestep mapping for the ids
|
||||
* @param[in] id_max The globally largest id used previously
|
||||
* @param[out] new_ids The newly renumbered blob ids (0:ids.max())
|
||||
*/
|
||||
|
|
@ -155,9 +159,9 @@ void getNewIDs( ID_map_struct& map, BlobIDType& id_max, std::vector<BlobIDType>&
|
|||
* @details This functions computes the map of blob ids between iterations.
|
||||
* Note: we also update the map to reflect the new ids
|
||||
* @param[out] new_ids The newly renumbered blob ids (0:ids.max())
|
||||
* @param[in/out] IDs The blob ids to renumber
|
||||
* @param[in,out] IDs The blob ids to renumber
|
||||
*/
|
||||
void renumberIDs( const std::vector<BlobIDType>& new_id_list, BlobIDArray& IDs );
|
||||
void renumberIDs( const std::vector<BlobIDType>& new_ids, BlobIDArray& IDs );
|
||||
|
||||
|
||||
/*!
|
||||
|
|
|
|||
|
|
@ -40,6 +40,7 @@ inline bool operator<(const Vec& l, const Vec& r){ return l.x*l.x+l.y*l.y+l.z*l.
|
|||
* @param[in] ID Segmentation id
|
||||
* @param[in] Dm Domain information
|
||||
* @param[in] periodic Directions that are periodic
|
||||
* @param[in] dx Cell size
|
||||
*/
|
||||
template<class TYPE>
|
||||
void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
|
||||
|
|
@ -52,6 +53,7 @@ void CalcDist( Array<TYPE> &Distance, const Array<char> &ID, const Domain &Dm,
|
|||
* @param[in] ID Domain id
|
||||
* @param[in] Dm Domain information
|
||||
* @param[in] periodic Directions that are periodic
|
||||
* @param[in] dx Cell size
|
||||
*/
|
||||
void CalcVecDist( Array<Vec> &Distance, const Array<int> &ID, const Domain &Dm,
|
||||
const std::array<bool,3>& periodic = {true,true,true}, const std::array<double,3>& dx = {1,1,1} );
|
||||
|
|
|
|||
|
|
@ -41,7 +41,10 @@ void Med3D( const Array<float> &Input, Array<float> &Output );
|
|||
* @details This routine performs a non-linear local means filter
|
||||
* @param[in] Input Input image
|
||||
* @param[in] Mean Mean value
|
||||
* @param[in] Distance Distance
|
||||
* @param[out] Output Output image
|
||||
* @param[in] d
|
||||
* @param[in] h
|
||||
*/
|
||||
int NLM3D( const Array<float> &Input, Array<float> &Mean,
|
||||
const Array<float> &Distance, Array<float> &Output, const int d, const float h);
|
||||
|
|
|
|||
|
|
@ -105,7 +105,7 @@ Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<Array<TY
|
|||
* multidimensional filter H. The result B has the same size and class as A.
|
||||
* This version works with separable filters and is more efficient than a single filter.
|
||||
* @param[in] A The input array (Nx,Ny,Nz)
|
||||
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
|
||||
* @param[in] Nh The filter size
|
||||
* @param[in] boundary The boundary conditions to apply (ndim):
|
||||
* fixed - Input array values outside the bounds of the array are
|
||||
* implicitly assumed to have the value X
|
||||
|
|
@ -130,6 +130,7 @@ Array<TYPE> imfilter_separable( const Array<TYPE>& A, const std::vector<int>& Nh
|
|||
* This version works with separable filters and is more efficient than a single filter.
|
||||
* @param[in] A The input array (Nx,Ny,Nz)
|
||||
* @param[in] H The filter [2*Nhx+1,2*Nhy+1,...]
|
||||
* @param[in] Nh The filter size
|
||||
* @param[in] boundary The boundary conditions to apply (ndim):
|
||||
* fixed - Input array values outside the bounds of the array are
|
||||
* implicitly assumed to have the value X
|
||||
|
|
|
|||
|
|
@ -67,6 +67,7 @@ class fillHalo
|
|||
public:
|
||||
/*!
|
||||
* @brief Default constructor
|
||||
* @param[in] comm Communicator to use
|
||||
* @param[in] info Rank and neighbor rank info
|
||||
* @param[in] n Number of local cells
|
||||
* @param[in] ng Number of ghost cells
|
||||
|
|
|
|||
|
|
@ -821,7 +821,7 @@ public: // Member functions
|
|||
* @return Output array for allGather
|
||||
*/
|
||||
template<class type>
|
||||
std::vector<type> allGather( const std::vector<type> &x_in ) const;
|
||||
std::vector<type> allGather( const std::vector<type> &x ) const;
|
||||
|
||||
|
||||
/*!
|
||||
|
|
|
|||
|
|
@ -148,6 +148,8 @@ extern "C" void ScaLBL_D3Q7_AAeven_Poisson_ElectricPotential(int *Map, double *d
|
|||
|
||||
extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, int start, int finish, int Np);
|
||||
|
||||
//extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish);
|
||||
|
||||
//maybe deprecated
|
||||
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
|
||||
// int strideY, int strideZ,int start, int finish, int Np);
|
||||
|
|
|
|||
|
|
@ -46,6 +46,7 @@ using StackTrace::Utilities::sleep_s;
|
|||
* \details This routine will peform the default startup sequence
|
||||
* \param argc argc from main
|
||||
* \param argv argv from main
|
||||
* \param multiple Intialize mpi with MPI_THREAD_MULTIPLE support?
|
||||
*/
|
||||
void startup( int argc, char **argv, bool multiple=true );
|
||||
|
||||
|
|
|
|||
|
|
@ -235,6 +235,87 @@ extern "C" void ScaLBL_D3Q7_Poisson_Init(int *Map, double *dist, double *Psi, in
|
|||
}
|
||||
}
|
||||
|
||||
extern "C" void ScaLBL_D3Q7_PoissonResidualError(int *neighborList, int *Map, double *ResidualError, double *Psi, double *Den_charge, double epsilon_LB,int strideY, int strideZ,int start, int finish){
|
||||
|
||||
int n,nn,ijk;
|
||||
double psi;//electric potential
|
||||
double rho_e;//local charge density
|
||||
// neighbors of electric potential psi
|
||||
double m1,m2,m4,m6,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18;
|
||||
double m0,m3,m5,m7;
|
||||
double psi_Laplacian;
|
||||
double residual_error;
|
||||
|
||||
for (n=start; n<finish; n++){
|
||||
|
||||
//Load data
|
||||
rho_e = Den_charge[n];
|
||||
ijk=Map[n];
|
||||
psi = Psi[ijk];
|
||||
|
||||
// COMPUTE THE COLOR GRADIENT
|
||||
//........................................................................
|
||||
//.................Read Phase Indicator Values............................
|
||||
//........................................................................
|
||||
nn = ijk-1; // neighbor index (get convention)
|
||||
m1 = Psi[nn]; // get neighbor for phi - 1
|
||||
//........................................................................
|
||||
nn = ijk+1; // neighbor index (get convention)
|
||||
m2 = Psi[nn]; // get neighbor for phi - 2
|
||||
//........................................................................
|
||||
nn = ijk-strideY; // neighbor index (get convention)
|
||||
m3 = Psi[nn]; // get neighbor for phi - 3
|
||||
//........................................................................
|
||||
nn = ijk+strideY; // neighbor index (get convention)
|
||||
m4 = Psi[nn]; // get neighbor for phi - 4
|
||||
//........................................................................
|
||||
nn = ijk-strideZ; // neighbor index (get convention)
|
||||
m5 = Psi[nn]; // get neighbor for phi - 5
|
||||
//........................................................................
|
||||
nn = ijk+strideZ; // neighbor index (get convention)
|
||||
m6 = Psi[nn]; // get neighbor for phi - 6
|
||||
//........................................................................
|
||||
nn = ijk-strideY-1; // neighbor index (get convention)
|
||||
m7 = Psi[nn]; // get neighbor for phi - 7
|
||||
//........................................................................
|
||||
nn = ijk+strideY+1; // neighbor index (get convention)
|
||||
m8 = Psi[nn]; // get neighbor for phi - 8
|
||||
//........................................................................
|
||||
nn = ijk+strideY-1; // neighbor index (get convention)
|
||||
m9 = Psi[nn]; // get neighbor for phi - 9
|
||||
//........................................................................
|
||||
nn = ijk-strideY+1; // neighbor index (get convention)
|
||||
m10 = Psi[nn]; // get neighbor for phi - 10
|
||||
//........................................................................
|
||||
nn = ijk-strideZ-1; // neighbor index (get convention)
|
||||
m11 = Psi[nn]; // get neighbor for phi - 11
|
||||
//........................................................................
|
||||
nn = ijk+strideZ+1; // neighbor index (get convention)
|
||||
m12 = Psi[nn]; // get neighbor for phi - 12
|
||||
//........................................................................
|
||||
nn = ijk+strideZ-1; // neighbor index (get convention)
|
||||
m13 = Psi[nn]; // get neighbor for phi - 13
|
||||
//........................................................................
|
||||
nn = ijk-strideZ+1; // neighbor index (get convention)
|
||||
m14 = Psi[nn]; // get neighbor for phi - 14
|
||||
//........................................................................
|
||||
nn = ijk-strideZ-strideY; // neighbor index (get convention)
|
||||
m15 = Psi[nn]; // get neighbor for phi - 15
|
||||
//........................................................................
|
||||
nn = ijk+strideZ+strideY; // neighbor index (get convention)
|
||||
m16 = Psi[nn]; // get neighbor for phi - 16
|
||||
//........................................................................
|
||||
nn = ijk+strideZ-strideY; // neighbor index (get convention)
|
||||
m17 = Psi[nn]; // get neighbor for phi - 17
|
||||
//........................................................................
|
||||
nn = ijk-strideZ+strideY; // neighbor index (get convention)
|
||||
m18 = Psi[nn]; // get neighbor for phi - 18
|
||||
|
||||
psi_Laplacian = 2.0*3.0/18.0*(m1+m2+m3+m4+m5+m6-6*psi+0.5*(m7+m8+m9+m10+m11+m12+m13+m14+m15+m16+m17+m18-12*psi));//Laplacian of electric potential
|
||||
residual_error = psi_Laplacian+rho_e/epsilon_LB;
|
||||
ResidualError[n] = residual_error;
|
||||
}
|
||||
}
|
||||
//extern "C" void ScaLBL_D3Q7_Poisson_ElectricField(int *neighborList, int *Map, signed char *ID, double *Psi, double *ElectricField, int SolidBC,
|
||||
// int strideY, int strideZ,int start, int finish, int Np){
|
||||
//
|
||||
|
|
|
|||
|
|
@ -379,12 +379,6 @@ MAX_INITIALIZER_LINES = 30
|
|||
|
||||
SHOW_USED_FILES = YES
|
||||
|
||||
# If the sources in your project are distributed over multiple directories
|
||||
# then setting the SHOW_DIRECTORIES tag to YES will show the directory hierarchy
|
||||
# in the documentation. The default is NO.
|
||||
|
||||
SHOW_DIRECTORIES = YES
|
||||
|
||||
# Set the SHOW_FILES tag to NO to disable the generation of the Files page.
|
||||
# This will remove the Files entry from the Quick Index and from the
|
||||
# Folder Tree View (if specified). The default is YES.
|
||||
|
|
@ -467,7 +461,7 @@ INPUT_ENCODING = UTF-8
|
|||
# *.hxx *.hpp *.h++ *.idl *.odl *.cs *.php *.php3 *.inc *.m *.mm *.dox *.py
|
||||
# *.f90 *.f *.for *.vhd *.vhdl
|
||||
|
||||
FILE_PATTERNS = *.h *.hh *.cu
|
||||
FILE_PATTERNS = *.h *.hh
|
||||
|
||||
# The RECURSIVE tag can be used to turn specify whether or not subdirectories
|
||||
# should be searched for input files as well. Possible values are YES and NO.
|
||||
|
|
@ -728,12 +722,6 @@ HTML_COLORSTYLE_GAMMA = 80
|
|||
|
||||
HTML_TIMESTAMP = YES
|
||||
|
||||
# If the HTML_ALIGN_MEMBERS tag is set to YES, the members of classes,
|
||||
# files or namespaces will be aligned in HTML using tables. If set to
|
||||
# NO a bullet list will be used.
|
||||
|
||||
HTML_ALIGN_MEMBERS = YES
|
||||
|
||||
# If the HTML_DYNAMIC_SECTIONS tag is set to YES then the generated HTML
|
||||
# documentation will contain sections that can be hidden and shown after the
|
||||
# page has loaded. For this to work a browser that supports
|
||||
|
|
@ -850,7 +838,7 @@ TREEVIEW_WIDTH = 250
|
|||
# If the GENERATE_LATEX tag is set to YES (the default) Doxygen will
|
||||
# generate Latex output.
|
||||
|
||||
GENERATE_LATEX = YES
|
||||
GENERATE_LATEX = NO
|
||||
|
||||
# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put.
|
||||
# If a relative path is entered the value of OUTPUT_DIRECTORY will be
|
||||
|
|
@ -1005,18 +993,6 @@ GENERATE_XML = NO
|
|||
|
||||
XML_OUTPUT = xml
|
||||
|
||||
# The XML_SCHEMA tag can be used to specify an XML schema,
|
||||
# which can be used by a validating XML parser to check the
|
||||
# syntax of the XML files.
|
||||
|
||||
XML_SCHEMA =
|
||||
|
||||
# The XML_DTD tag can be used to specify an XML DTD,
|
||||
# which can be used by a validating XML parser to check the
|
||||
# syntax of the XML files.
|
||||
|
||||
XML_DTD =
|
||||
|
||||
# If the XML_PROGRAMLISTING tag is set to YES Doxygen will
|
||||
# dump the program listings (including syntax highlighting
|
||||
# and cross-referencing information) to the XML output. Note that
|
||||
|
|
@ -1174,11 +1150,6 @@ ALLEXTERNALS = NO
|
|||
|
||||
EXTERNAL_GROUPS = YES
|
||||
|
||||
# The PERL_PATH should be the absolute path and name of the perl script
|
||||
# interpreter (i.e. the result of `which perl').
|
||||
|
||||
PERL_PATH = /usr/bin/perl
|
||||
|
||||
#---------------------------------------------------------------------------
|
||||
# Configuration options related to the dot tool
|
||||
#---------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -4,8 +4,6 @@
|
|||
*
|
||||
* - \ref IO "IO routines"
|
||||
* - \ref Utilities "Utility routines"
|
||||
* - \ref silo "Access to silo routines"
|
||||
* - \ref netcdf "Access to netcdf routines"
|
||||
*
|
||||
* \author James McClure
|
||||
*/
|
||||
|
|
|
|||
|
|
@ -627,13 +627,13 @@ void ScaLBL_IonModel::AssignSolidBoundary(double *ion_solid)
|
|||
}
|
||||
}
|
||||
|
||||
void ScaLBL_IonModel::AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion)
|
||||
void ScaLBL_IonModel::AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion, int ic)
|
||||
{
|
||||
double *Ci_host;
|
||||
Ci_host = new double[N];
|
||||
double VALUE=0.f;
|
||||
|
||||
Mask->ReadFromFile(File_ion[0],File_ion[1],Ci_host);
|
||||
Mask->ReadFromFile(File_ion[2*ic],File_ion[2*ic+1],Ci_host);
|
||||
|
||||
for (int k=0;k<Nz;k++){
|
||||
for (int j=0;j<Ny;j++){
|
||||
|
|
@ -722,20 +722,24 @@ void ScaLBL_IonModel::Initialize(){
|
|||
*/
|
||||
if (rank==0) printf ("LB Ion Solver: initializing D3Q7 distributions\n");
|
||||
if (ion_db->keyExists("IonConcentrationFile")){
|
||||
//TODO: Need to figure out how to deal with multi-species concentration initialization
|
||||
//NOTE: "IonConcentrationFile" is a vector, including "file_name, datatype"
|
||||
auto File_ion = ion_db->getVector<std::string>( "IonConcentrationFile" );
|
||||
double *Ci_host;
|
||||
Ci_host = new double[number_ion_species*Np];
|
||||
for (size_t ic=0; ic<number_ion_species; ic++){
|
||||
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion);
|
||||
if (File_ion.size()==2*number_ion_species){
|
||||
double *Ci_host;
|
||||
Ci_host = new double[number_ion_species*Np];
|
||||
for (size_t ic=0; ic<number_ion_species; ic++){
|
||||
AssignIonConcentration_FromFile(&Ci_host[ic*Np],File_ion,ic);
|
||||
}
|
||||
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
|
||||
comm.barrier();
|
||||
for (size_t ic=0; ic<number_ion_species; ic++){
|
||||
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
|
||||
}
|
||||
delete [] Ci_host;
|
||||
}
|
||||
ScaLBL_CopyToDevice(Ci, Ci_host, number_ion_species*sizeof(double)*Np);
|
||||
comm.barrier();
|
||||
for (size_t ic=0; ic<number_ion_species; ic++){
|
||||
ScaLBL_D3Q7_Ion_Init_FromFile(&fq[ic*Np*7],&Ci[ic*Np],Np);
|
||||
else{
|
||||
ERROR("Error: Number of user-input ion concentration files should be equal to number of ion species!\n");
|
||||
}
|
||||
delete [] Ci_host;
|
||||
}
|
||||
else{
|
||||
for (size_t ic=0; ic<number_ion_species; ic++){
|
||||
|
|
|
|||
|
|
@ -96,7 +96,7 @@ private:
|
|||
//int rank,nprocs;
|
||||
void LoadParams(std::shared_ptr<Database> db0);
|
||||
void AssignSolidBoundary(double *ion_solid);
|
||||
void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion);
|
||||
void AssignIonConcentration_FromFile(double *Ci,const vector<std::string> &File_ion,int ic);
|
||||
void IonConcentration_LB_to_Phys(DoubleArray &Den_reg);
|
||||
};
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -69,6 +69,8 @@ void ScaLBL_Poisson::ReadParams(string filename){
|
|||
if (electric_db->keyExists( "tolerance" )){
|
||||
tolerance = electric_db->getScalar<double>( "tolerance" );
|
||||
}
|
||||
//'tolerance_method' can be {"MSE","MSE_max"}
|
||||
tolerance_method = electric_db->getWithDefault<std::string>( "tolerance_method", "MSE" );
|
||||
if (electric_db->keyExists( "epsilonR" )){
|
||||
epsilonR = electric_db->getScalar<double>( "epsilonR" );
|
||||
}
|
||||
|
|
@ -122,6 +124,15 @@ void ScaLBL_Poisson::ReadParams(string filename){
|
|||
if (rank==0) printf("LB-Poisson Solver: steady-state MaxTimeStep = %i; steady-state tolerance = %.3g \n", timestepMax,tolerance);
|
||||
if (rank==0) printf(" LB relaxation tau = %.5g \n", tau);
|
||||
if (rank==0) printf("***********************************************************************************\n");
|
||||
if (tolerance_method.compare("MSE")==0){
|
||||
if (rank==0) printf("LB-Poisson Solver: Use averaged MSE to check solution convergence.\n");
|
||||
}
|
||||
else if (tolerance_method.compare("MSE_max")==0){
|
||||
if (rank==0) printf("LB-Poisson Solver: Use maximum MSE to check solution convergence.\n");
|
||||
}
|
||||
else{
|
||||
if (rank==0) printf("LB-Poisson Solver: tolerance_method=%s cannot be identified!\n",tolerance_method.c_str());
|
||||
}
|
||||
|
||||
switch (BoundaryConditionSolid){
|
||||
case 1:
|
||||
|
|
@ -150,6 +161,7 @@ void ScaLBL_Poisson::SetDomain(){
|
|||
N = Nx*Ny*Nz;
|
||||
Distance.resize(Nx,Ny,Nz);
|
||||
Psi_host.resize(Nx,Ny,Nz);
|
||||
Psi_previous.resize(Nx,Ny,Nz);
|
||||
|
||||
for (int i=0; i<Nx*Ny*Nz; i++) Dm->id[i] = 1; // initialize this way
|
||||
//Averages = std::shared_ptr<TwoPhase> ( new TwoPhase(Dm) ); // TwoPhase analysis object
|
||||
|
|
@ -338,6 +350,7 @@ void ScaLBL_Poisson::Create(){
|
|||
ScaLBL_AllocateDeviceMemory((void **) &fq, 7*dist_mem_size);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &Psi, sizeof(double)*Nx*Ny*Nz);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &ElectricField, 3*sizeof(double)*Np);
|
||||
ScaLBL_AllocateDeviceMemory((void **) &ResidualError, sizeof(double)*Np);
|
||||
//...........................................................................
|
||||
|
||||
// Update GPU data structures
|
||||
|
|
@ -541,12 +554,10 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
|
|||
|
||||
timestep=0;
|
||||
double error = 1.0;
|
||||
double psi_avg_previous = 0.0;
|
||||
while (timestep < timestepMax && error > tolerance) {
|
||||
//************************************************************************/
|
||||
// *************ODD TIMESTEP*************//
|
||||
timestep++;
|
||||
|
||||
SolveElectricPotentialAAodd(timestep_from_Study);//update electric potential
|
||||
SolvePoissonAAodd(ChargeDensity);//perform collision
|
||||
ScaLBL_Comm->Barrier(); comm.barrier();
|
||||
|
|
@ -559,33 +570,86 @@ void ScaLBL_Poisson::Run(double *ChargeDensity, int timestep_from_Study){
|
|||
//************************************************************************/
|
||||
|
||||
// Check convergence of steady-state solution
|
||||
if (timestep==2){
|
||||
//save electric potential for convergence check
|
||||
ScaLBL_CopyToHost(Psi_previous.data(),Psi,sizeof(double)*Nx*Ny*Nz);
|
||||
}
|
||||
if (timestep%analysis_interval==0){
|
||||
|
||||
//ScaLBL_Comm->RegularLayout(Map,Psi,Psi_host);
|
||||
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
|
||||
double count_loc=0;
|
||||
double count;
|
||||
double psi_avg;
|
||||
double psi_loc=0.f;
|
||||
if (tolerance_method.compare("MSE")==0){
|
||||
double count_loc=0;
|
||||
double count;
|
||||
double MSE_loc=0.0;
|
||||
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
if (Distance(i,j,k) > 0){
|
||||
MSE_loc += (Psi_host(i,j,k) - Psi_previous(i,j,k))*(Psi_host(i,j,k) - Psi_previous(i,j,k));
|
||||
count_loc+=1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
error=Dm->Comm.sumReduce(MSE_loc);
|
||||
count=Dm->Comm.sumReduce(count_loc);
|
||||
error /= count;
|
||||
}
|
||||
else if (tolerance_method.compare("MSE_max")==0){
|
||||
vector<double>MSE_loc;
|
||||
double MSE_loc_max;
|
||||
ScaLBL_CopyToHost(Psi_host.data(),Psi,sizeof(double)*Nx*Ny*Nz);
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
if (Distance(i,j,k) > 0){
|
||||
MSE_loc.push_back((Psi_host(i,j,k) - Psi_previous(i,j,k))*(Psi_host(i,j,k) - Psi_previous(i,j,k)));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
vector<double>::iterator it_max = max_element(MSE_loc.begin(),MSE_loc.end());
|
||||
unsigned int idx_max=distance(MSE_loc.begin(),it_max);
|
||||
MSE_loc_max=MSE_loc[idx_max];
|
||||
error=Dm->Comm.maxReduce(MSE_loc_max);
|
||||
}
|
||||
else{
|
||||
ERROR("Error: user-specified tolerance_method cannot be identified; check you input database! \n");
|
||||
}
|
||||
ScaLBL_CopyToHost(Psi_previous.data(),Psi,sizeof(double)*Nx*Ny*Nz);
|
||||
|
||||
for (int k=1; k<Nz-1; k++){
|
||||
for (int j=1; j<Ny-1; j++){
|
||||
for (int i=1; i<Nx-1; i++){
|
||||
if (Distance(i,j,k) > 0){
|
||||
psi_loc += Psi_host(i,j,k);
|
||||
count_loc+=1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
psi_avg=Dm->Comm.sumReduce( psi_loc);
|
||||
count=Dm->Comm.sumReduce( count_loc);
|
||||
|
||||
psi_avg /= count;
|
||||
double psi_avg_mag=psi_avg;
|
||||
if (psi_avg==0.0) psi_avg_mag=1.0;
|
||||
error = fabs(psi_avg-psi_avg_previous)/fabs(psi_avg_mag);
|
||||
psi_avg_previous = psi_avg;
|
||||
|
||||
|
||||
|
||||
//legacy code that tried to use residual to check convergence
|
||||
//ScaLBL_D3Q7_PoissonResidualError(NeighborList,dvcMap,ResidualError,Psi,ChargeDensity,epsilon_LB,Nx,Nx*Ny,ScaLBL_Comm->FirstInterior(), ScaLBL_Comm->LastInterior());
|
||||
//ScaLBL_D3Q7_PoissonResidualError(NeighborList,dvcMap,ResidualError,Psi,ChargeDensity,epsilon_LB,Nx,Nx*Ny,0, ScaLBL_Comm->LastExterior());
|
||||
//ScaLBL_Comm->Barrier(); comm.barrier();
|
||||
|
||||
//vector<double> ResidualError_host(Np);
|
||||
//double error_loc_max;
|
||||
////calculate the maximum residual error
|
||||
//ScaLBL_CopyToHost(&ResidualError_host[0],ResidualError,sizeof(double)*Np);
|
||||
|
||||
//vector<double>::iterator it_temp1,it_temp2;
|
||||
//it_temp1=ResidualError_host.begin();
|
||||
//advance(it_temp1,ScaLBL_Comm->LastExterior());
|
||||
//vector<double>::iterator it_max = max_element(ResidualError_host.begin(),it_temp1);
|
||||
//unsigned int idx_max1 = distance(ResidualError_host.begin(),it_max);
|
||||
|
||||
//it_temp1=ResidualError_host.begin();
|
||||
//it_temp2=ResidualError_host.begin();
|
||||
//advance(it_temp1,ScaLBL_Comm->FirstInterior());
|
||||
//advance(it_temp2,ScaLBL_Comm->LastInterior());
|
||||
//it_max = max_element(it_temp1,it_temp2);
|
||||
//unsigned int idx_max2 = distance(ResidualError_host.begin(),it_max);
|
||||
//if (ResidualError_host[idx_max1]>ResidualError_host[idx_max2]){
|
||||
// error_loc_max=ResidualError_host[idx_max1];
|
||||
//}
|
||||
//else{
|
||||
// error_loc_max=ResidualError_host[idx_max2];
|
||||
//}
|
||||
//error = Dm->Comm.maxReduce(error_loc_max);
|
||||
}
|
||||
}
|
||||
if(WriteLog==true){
|
||||
|
|
|
|||
|
|
@ -10,6 +10,7 @@
|
|||
#include <stdexcept>
|
||||
#include <fstream>
|
||||
#include <cmath>
|
||||
#include <iterator>
|
||||
|
||||
#include "common/ScaLBL.h"
|
||||
#include "common/Communication.h"
|
||||
|
|
@ -48,6 +49,7 @@ public:
|
|||
int BoundaryConditionSolid;
|
||||
double tau;
|
||||
double tolerance;
|
||||
std::string tolerance_method;
|
||||
double k2_inv;
|
||||
double epsilon0,epsilon0_LB,epsilonR,epsilon_LB;
|
||||
double Vin, Vout;
|
||||
|
|
@ -78,6 +80,7 @@ public:
|
|||
IntArray Map;
|
||||
DoubleArray Distance;
|
||||
DoubleArray Psi_host;
|
||||
DoubleArray Psi_previous;
|
||||
int *NeighborList;
|
||||
int *dvcMap;
|
||||
//signed char *dvcID;
|
||||
|
|
@ -85,6 +88,7 @@ public:
|
|||
double *Psi;
|
||||
double *ElectricField;
|
||||
double *ChargeDensityDummy;// for debugging
|
||||
double *ResidualError;
|
||||
|
||||
private:
|
||||
Utilities::MPI comm;
|
||||
|
|
|
|||
|
|
@ -237,7 +237,6 @@ inline int32_atomic atomic_fetch_and_or( int32_atomic volatile *v, int32_atomic
|
|||
* \brief Fetch the current value and "ou" with given value
|
||||
* \details Perform *v = (*v) | x, returning the previous value
|
||||
* \return Returns the previous value before the "and" operation
|
||||
* \param[in] v The pointer to the value to check and swap
|
||||
* \param[in] v The pointer to the value to check and or
|
||||
* \param[in] x The value to or
|
||||
*/
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user