1#ifndef RPC_LR_COMPRESSOR_TPP
2#define RPC_LR_COMPRESSOR_TPP
11#include "LrCompressor.h"
12#include <rpc/System.h>
13#include <rpc/fts/compressor/intra/IntraCorrelation.h>
15#include <util/format/Dbl.h>
16#include <pscf/mesh/MeshIterator.h>
17#include <pscf/iterator/NanException.h>
18#include <prdc/crystal/shiftToMinimum.h>
34 errorType_(
"rmsResid"),
37 isIntraCalculated_(false),
51 read(in,
"epsilon", epsilon_);
52 readOptional(in,
"maxItr", maxItr_);
53 readOptional(in,
"verbose", verbose_);
54 readOptional(in,
"errorType", errorType_);
61 const int nMonomer = system().mixture().nMonomer();
62 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
63 for (
int i = 0; i < D; ++i) {
65 kMeshDimensions_[i] = dimensions[i];
67 kMeshDimensions_[i] = dimensions[i]/2 + 1;
73 resid_.allocate(dimensions);
74 residK_.allocate(dimensions);
75 wFieldTmp_.allocate(nMonomer);
76 intraCorrelationK_.allocate(kMeshDimensions_);
77 for (
int i = 0; i < nMonomer; ++i) {
78 wFieldTmp_[i].allocate(dimensions);
84 if (!isIntraCalculated_){
85 intra_.computeIntraCorrelations(intraCorrelationK_);
86 isIntraCalculated_ =
true;
107 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
110 Log::file() <<
"------------------------------- \n";
121 error = computeError(verbose_);
123 Log::file() <<
", error = NaN" << std::endl;
127 Log::file() <<
", error = " <<
Dbl(error, 15) << std::endl;
131 if (error < epsilon_) {
137 Log::file() <<
"-------------------------------\n";
167 Log::file() <<
"Iterator failed to converge.\n";
176 const int nMonomer = system().mixture().nMonomer();
177 const int meshSize = system().domain().mesh().size();
180 for (
int i = 0 ; i < meshSize; ++i) {
185 for (
int j = 0; j < nMonomer; ++j) {
186 for (
int k = 0; k < meshSize; ++k) {
187 resid_[k] += system().c().rgrid(j)[k];
194 void LrCompressor<D>::updateWFields()
196 const int nMonomer = system().mixture().nMonomer();
197 const int meshSize = system().domain().mesh().size();
198 const double vMonomer = system().mixture().vMonomer();
201 system().domain().fft().forwardTransform(resid_, residK_);
202 MeshIterator<D> iter;
203 iter.setDimensions(residK_.dftDimensions());
204 for (iter.begin(); !iter.atEnd(); ++iter) {
205 residK_[iter.rank()][0] *= 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
206 residK_[iter.rank()][1] *= 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
210 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
212 for (
int i = 0; i < nMonomer; i++){
213 for (
int k = 0; k < meshSize; k++){
214 wFieldTmp_[i][k] = system().w().rgrid(i)[k] + resid_[k];
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";
256 for (
int i = 0; i < n; i++) {
258 if (std::isnan(value)) {
259 throw NanException(
"LrCompressor::dotProduct",__FILE__,__LINE__,0);
261 if (fabs(value) > max) {
270 double LrCompressor<D>::dotProduct(RField<D>
const & a,
273 const int n = a.capacity();
276 for (
int i = 0; i < n; i++) {
278 if (std::isnan(a[i]) || std::isnan(b[i])) {
279 throw NanException(
"AmCompressor::dotProduct",__FILE__,__LINE__,0);
288 double LrCompressor<D>::norm(RField<D>
const & a)
290 double normSq = dotProduct(a, a);
296 double LrCompressor<D>::computeError(
int verbose)
298 const int meshSize = system().domain().mesh().size();
302 double maxRes =
maxAbs(resid_);
304 double normRes = norm(resid_);
306 double rmsRes = normRes/sqrt(meshSize);
307 if (errorType_ ==
"maxResid") {
309 }
else if (errorType_ ==
"normResid") {
311 }
else if (errorType_ ==
"rmsResid") {
314 UTIL_THROW(
"Invalid iterator error type in parameter file.");
319 Log::file() <<
"Max Residual = " <<
Dbl(maxRes,15) <<
"\n";
320 Log::file() <<
"Residual Norm = " <<
Dbl(normRes,15) <<
"\n";
321 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 setClassName(const char *className)
Set class name string.
LrCompressor(System< D > &system)
Constructor.
~LrCompressor()
Destructor.
int compress()
Iterate pressure field to obtain partial saddle point.
void outputTimers(std::ostream &out)
Return compressor times contributions.
void clearTimers()
Clear timers.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setup()
Initialize just before entry to iterative loop.
Main class for SCFT or PS-FTS simulation of one system.
int capacity() const
Return allocated size.
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 maxAbs(DeviceArray< cudaReal > const &in)
Get maximum absolute magnitude of array elements (GPU kernel wrapper).
cudaReal max(DeviceArray< cudaReal > const &in)
Get maximum of array elements (GPU kernel wrapper).
PSCF package top-level namespace.
Utility classes for scientific computation.
float product(float a, float b)
Product for float Data.