1#ifndef RPG_LR_COMPRESSOR_TPP
2#define RPG_LR_COMPRESSOR_TPP
11#include "LrCompressor.h"
12#include <rpg/solvers/Mixture.h>
13#include <rpg/field/Domain.h>
14#include <rpg/system/System.h>
15#include <prdc/crystal/shiftToMinimum.h>
16#include <prdc/cuda/resources.h>
17#include <pscf/mesh/MeshIterator.h>
18#include <pscf/iterator/NanException.h>
19#include <util/format/Dbl.h>
40 isIntraCalculated_(false),
63 read(in,
"epsilon", epsilon_);
75 const int nMonomer =
system().mixture().nMonomer();
82 resid_.allocate(dimensions);
83 residK_.allocate(dimensions);
84 wFieldTmp_.allocate(nMonomer);
85 intraCorrelationK_.allocate(kMeshDimensions_);
86 for (
int i = 0; i < nMonomer; ++i) {
87 wFieldTmp_[i].allocate(dimensions);
93 if (!isIntraCalculated_){
94 intra_.computeIntraCorrelations(intraCorrelationK_);
95 isIntraCalculated_ =
true;
120 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
123 Log::file() <<
"------------------------------- \n";
134 error = computeError(verbose_);
136 Log::file() <<
", error = NaN" << std::endl;
140 Log::file() <<
", error = " <<
Dbl(error, 15) << std::endl;
144 if (error < epsilon_) {
150 Log::file() <<
"-------------------------------\n";
180 Log::file() <<
"Iterator failed to converge.\n";
189 void LrCompressor<D>::getResidual()
191 const int nMonomer = system().mixture().nMonomer();
197 for (
int i = 1; i < nMonomer; i++) {
206 void LrCompressor<D>::updateWFields(){
207 const int nMonomer = system().mixture().nMonomer();
208 const double vMonomer = system().mixture().vMonomer();
211 system().domain().fft().forwardTransform(resid_, residK_);
217 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
220 for (
int i = 0; i < nMonomer; i++) {
221 VecOp::addVV(wFieldTmp_[i], system().w().rgrid(i), resid_);
225 system().w().setRGrid(wFieldTmp_);
232 void LrCompressor<D>::outputToLog()
242 double total = timerTotal_.time();
244 out <<
"LrCompressor time contributions:\n";
246 out <<
"Total" << std::setw(22)<<
"Per Iteration"
247 << std::setw(9) <<
"Fraction" <<
"\n";
248 out <<
"MDE solution: "
249 <<
Dbl(timerMDE_.time(), 9, 3) <<
" s, "
250 <<
Dbl(timerMDE_.time()/totalItr_, 9, 3) <<
" s, "
251 <<
Dbl(timerMDE_.time()/total, 9, 3) <<
"\n";
268 double LrCompressor<D>::maxAbs(
RField<D> const & a)
277 double LrCompressor<D>::dotProduct(RField<D>
const & a,
288 double LrCompressor<D>::norm(RField<D>
const & a)
290 double normSq = dotProduct(a, a);
298 double LrCompressor<D>::computeError(
int verbose)
300 const int meshSize = system().domain().mesh().size();
304 double maxRes =
maxAbs(resid_);
306 double normRes = norm(resid_);
308 double rmsRes = normRes/sqrt(meshSize);
309 if (errorType_ ==
"max") {
311 }
else if (errorType_ ==
"norm") {
313 }
else if (errorType_ ==
"rms") {
316 UTIL_THROW(
"Invalid iterator error type in parameter file.");
321 Log::file() <<
"Max Residual = " <<
Dbl(maxRes,15) <<
"\n";
322 Log::file() <<
"Residual Norm = " <<
Dbl(normRes,15) <<
"\n";
323 Log::file() <<
"RMS Residual = " <<
Dbl(rmsRes,15) <<
"\n";
An IntVec<D, T> is a D-component vector of elements of integer type T.
Exception thrown when not-a-number (NaN) is encountered.
static void computeKMesh(IntVec< D > const &rMeshDimensions, IntVec< D > &kMeshDimensions, int &kSize)
Compute dimensions and size of k-space mesh for DFT of real data.
Field of real double precision values on an FFT mesh.
int mdeCounter_
Count how many times MDE has been solved.
Compressor()
Default constructor.
System< D > const & system() const
Return const reference to parent system.
void setup()
Initialize just before entry to iterative loop.
void readParameters(std::istream &in)
Read all parameters and initialize.
ScalarParam< Type > & read(std::istream &in, const char *label, Type &value)
Add and read a new required ScalarParam < Type > object.
void setClassName(const char *className)
Set class name string.
void clearTimers()
Clear timers.
ScalarParam< Type > & readOptional(std::istream &in, const char *label, Type &value)
Add and read a new optional ScalarParam < Type > object.
~LrCompressor()
Destructor.
int compress()
Iterate pressure field to obtain partial saddle point.
LrCompressor(System< D > &system)
Constructor.
void outputTimers(std::ostream &out) const
Return compressor times contributions.
Main class, representing a complete physical system.
Wrapper for a double precision number, for formatted ostream output.
Wrapper for an int, for formatted ostream output.
static std::ostream & file()
Get log ostream by reference.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
#define UTIL_THROW(msg)
Macro for throwing an Exception, reporting function, file and line number.
double maxAbs(Array< double > const &in)
Get maximum absolute magnitude of array elements .
double innerProduct(Array< double > const &a, Array< double > const &b)
Compute inner product of two real arrays .
void addEqV(Array< double > &a, Array< double > const &b)
Vector addition in-place, a[i] += b[i].
void addVV(Array< double > &a, Array< double > const &b, Array< double > const &c)
Vector addition, a[i] = b[i] + c[i].
void subVS(Array< double > &a, Array< double > const &b, double c)
Vector subtraction, a[i] = b[i] - c.
void divEqVc(DeviceArray< cudaComplex > &a, DeviceArray< cudaReal > const &b, cudaReal const c)
Vector division in-place w/ coeff., a[i] /= (b[i] * c), kernel wrapper.
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.