1#ifndef RPG_LR_COMPRESSOR_TPP
2#define RPG_LR_COMPRESSOR_TPP
11#include "LrCompressor.h"
12#include <rpg/System.h>
13#include <rpg/fts/compressor/intra/IntraCorrelation.h>
14#include <prdc/crystal/shiftToMinimum.h>
15#include <prdc/cuda/resources.h>
16#include <pscf/mesh/MeshIterator.h>
17#include <pscf/iterator/NanException.h>
19#include <util/format/Dbl.h>
35 errorType_(
"rmsResid"),
38 isIntraCalculated_(false),
53 read(in,
"epsilon", epsilon_);
54 readOptional(in,
"maxItr", maxItr_);
55 readOptional(in,
"verbose", verbose_);
56 readOptional(in,
"errorType", errorType_);
63 const int nMonomer = system().mixture().nMonomer();
64 IntVec<D> const & dimensions = system().mesh().dimensions();
65 for (
int i = 0; i < D; ++i) {
67 kMeshDimensions_[i] = dimensions[i];
69 kMeshDimensions_[i] = dimensions[i]/2 + 1;
75 for (
int i = 0; i < D; ++i) {
76 kSize_ *= kMeshDimensions_[i];
81 resid_.allocate(dimensions);
82 residK_.allocate(dimensions);
83 wFieldTmp_.allocate(nMonomer);
84 intraCorrelationK_.allocate(kMeshDimensions_);
85 for (
int i = 0; i < nMonomer; ++i) {
86 wFieldTmp_[i].allocate(dimensions);
92 if (!isIntraCalculated_){
93 intra_.computeIntraCorrelations(intraCorrelationK_);
94 isIntraCalculated_ =
true;
116 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
119 Log::file() <<
"------------------------------- \n";
130 error = computeError(verbose_);
132 Log::file() <<
", error = NaN" << std::endl;
136 Log::file() <<
", error = " <<
Dbl(error, 15) << std::endl;
140 if (error < epsilon_) {
146 Log::file() <<
"-------------------------------\n";
176 Log::file() <<
"Iterator failed to converge.\n";
185 const int nMonomer = system().mixture().nMonomer();
191 for (
int i = 1; i < nMonomer; i++) {
198 void LrCompressor<D>::updateWFields(){
199 const int nMonomer = system().mixture().nMonomer();
200 const double vMonomer = system().mixture().vMonomer();
203 system().fft().forwardTransform(resid_, residK_);
209 system().fft().inverseTransformUnsafe(residK_, resid_);
212 for (
int i = 0; i < nMonomer; i++) {
213 VecOp::addVV(wFieldTmp_[i], system().w().rgrid(i), resid_);
217 system().setWRGrid(wFieldTmp_);
221 void LrCompressor<D>::outputToLog()
228 double total = timerTotal_.time();
230 out <<
"LrCompressor time contributions:\n";
232 out <<
"Total" << std::setw(22)<<
"Per Iteration"
233 << std::setw(9) <<
"Fraction" <<
"\n";
234 out <<
"MDE solution: "
235 <<
Dbl(timerMDE_.time(), 9, 3) <<
" s, "
236 <<
Dbl(timerMDE_.time()/totalItr_, 9, 3) <<
" s, "
237 <<
Dbl(timerMDE_.time()/total, 9, 3) <<
"\n";
258 double LrCompressor<D>::dotProduct(RField<D>
const & a,
267 double LrCompressor<D>::norm(RField<D>
const & a)
269 double normSq = dotProduct(a, a);
275 double LrCompressor<D>::computeError(
int verbose)
277 const int meshSize = system().domain().mesh().size();
281 double maxRes =
maxAbs(resid_);
283 double normRes = norm(resid_);
285 double rmsRes = normRes/sqrt(meshSize);
286 if (errorType_ ==
"maxResid") {
288 }
else if (errorType_ ==
"normResid") {
290 }
else if (errorType_ ==
"rmsResid") {
293 UTIL_THROW(
"Invalid iterator error type in parameter file.");
298 Log::file() <<
"Max Residual = " <<
Dbl(maxRes,15) <<
"\n";
299 Log::file() <<
"Residual Norm = " <<
Dbl(normRes,15) <<
"\n";
300 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.
Field of real double precision values on an FFT mesh.
Base class for iterators that impose incompressibility.
Linear response compressor.
void setup()
Initialize just before entry to iterative loop.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setClassName(const char *className)
Set class name string.
void clearTimers()
Clear timers.
~LrCompressor()
Destructor.
int compress()
Iterate pressure field to obtain partial saddle point.
void outputTimers(std::ostream &out)
Return compressor times contributions.
LrCompressor(System< D > &system)
Constructor.
Main class for calculations that represent one 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.
cudaReal innerProduct(DeviceArray< cudaReal > const &a, DeviceArray< cudaReal > const &b)
Compute inner product of two real arrays (GPU kernel wrapper).
cudaReal maxAbs(DeviceArray< cudaReal > const &in)
Get maximum absolute magnitude of array elements (GPU kernel wrapper).
void addEqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector addition in-place, a[i] += b[i], kernel wrapper (cudaReal).
void subVS(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, const int beginIdA, const int beginIdB, const int n)
Vector subtraction, a[i] = b[i] - c, kernel wrapper (cudaReal).
void addVV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, DeviceArray< cudaReal > const &c, const int beginIdA, const int beginIdB, const int beginIdC, const int n)
Vector addition, a[i] = b[i] + c[i], kernel wrapper (cudaReal).
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.
PSCF package top-level namespace.
Utility classes for scientific computation.