1#ifndef RPC_LR_COMPRESSOR_TPP
2#define RPC_LR_COMPRESSOR_TPP
11#include "LrCompressor.h"
12#include <rpc/system/System.h>
13#include <prdc/crystal/shiftToMinimum.h>
14#include <pscf/mesh/MeshIterator.h>
15#include <pscf/iterator/NanException.h>
17#include <util/format/Dbl.h>
33 errorType_(
"rmsResid"),
36 isIntraCalculated_(false),
50 read(in,
"epsilon", epsilon_);
60 const int nMonomer =
system().mixture().nMonomer();
62 for (
int i = 0; i < D; ++i) {
64 kMeshDimensions_[i] = dimensions[i];
66 kMeshDimensions_[i] = dimensions[i]/2 + 1;
72 resid_.allocate(dimensions);
73 residK_.allocate(dimensions);
74 wFieldTmp_.allocate(nMonomer);
75 intraCorrelationK_.allocate(kMeshDimensions_);
76 for (
int i = 0; i < nMonomer; ++i) {
77 wFieldTmp_[i].allocate(dimensions);
83 if (!isIntraCalculated_){
84 intra_.computeIntraCorrelations(intraCorrelationK_);
85 isIntraCalculated_ =
true;
109 for (itr_ = 0; itr_ < maxItr_; ++itr_) {
112 Log::file() <<
"------------------------------- \n";
123 error = computeError(verbose_);
125 Log::file() <<
", error = NaN" << std::endl;
129 Log::file() <<
", error = " <<
Dbl(error, 15) << std::endl;
133 if (error < epsilon_) {
139 Log::file() <<
"-------------------------------\n";
169 Log::file() <<
"Iterator failed to converge.\n";
178 void LrCompressor<D>::getResidual()
180 const int nMonomer = system().mixture().nMonomer();
181 const int meshSize = system().domain().mesh().size();
184 for (
int i = 0 ; i < meshSize; ++i) {
189 for (
int j = 0; j < nMonomer; ++j) {
190 for (
int k = 0; k < meshSize; ++k) {
191 resid_[k] += system().c().rgrid(j)[k];
198 void LrCompressor<D>::updateWFields()
200 const int nMonomer = system().mixture().nMonomer();
201 const int meshSize = system().domain().mesh().size();
202 const double vMonomer = system().mixture().vMonomer();
205 system().domain().fft().forwardTransform(resid_, residK_);
208 for (iter.begin(); !iter.atEnd(); ++iter) {
209 residK_[iter.rank()][0] *= 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
210 residK_[iter.rank()][1] *= 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
214 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
216 for (
int i = 0; i < nMonomer; i++){
217 for (
int k = 0; k < meshSize; k++){
218 wFieldTmp_[i][k] = system().w().rgrid(i)[k] + resid_[k];
221 system().w().setRGrid(wFieldTmp_);
225 void LrCompressor<D>::outputToLog()
232 double total = timerTotal_.time();
234 out <<
"LrCompressor time contributions:\n";
236 out <<
"Total" << std::setw(22)<<
"Per Iteration"
237 << std::setw(9) <<
"Fraction" <<
"\n";
238 out <<
"MDE solution: "
239 <<
Dbl(timerMDE_.time(), 9, 3) <<
" s, "
240 <<
Dbl(timerMDE_.time()/totalItr_, 9, 3) <<
" s, "
241 <<
Dbl(timerMDE_.time()/total, 9, 3) <<
"\n";
255 double LrCompressor<D>::maxAbs(
RField<D> const & a)
260 for (
int i = 0; i < n; i++) {
262 if (std::isnan(value)) {
263 throw NanException(
"LrCompressor::dotProduct",__FILE__,__LINE__,0);
265 if (fabs(value) > max) {
274 double LrCompressor<D>::dotProduct(
RField<D> const & a,
277 const int n = a.capacity();
280 for (
int i = 0; i < n; i++) {
282 if (std::isnan(a[i]) || std::isnan(b[i])) {
283 throw NanException(
"AmCompressor::dotProduct",__FILE__,__LINE__,0);
292 double LrCompressor<D>::norm(
RField<D> const & a)
294 double normSq = dotProduct(a, a);
300 double LrCompressor<D>::computeError(
int verbose)
302 const int meshSize = system().domain().mesh().size();
306 double maxRes = maxAbs(resid_);
308 double normRes = norm(resid_);
310 double rmsRes = normRes/sqrt(meshSize);
311 if (errorType_ ==
"maxResid") {
313 }
else if (errorType_ ==
"normResid") {
315 }
else if (errorType_ ==
"rmsResid") {
318 UTIL_THROW(
"Invalid iterator error type in parameter file.");
323 Log::file() <<
"Max Residual = " <<
Dbl(maxRes,15) <<
"\n";
324 Log::file() <<
"Residual Norm = " <<
Dbl(normRes,15) <<
"\n";
325 Log::file() <<
"RMS Residual = " <<
Dbl(rmsRes,15) <<
"\n";
An IntVec<D, T> is a D-component vector of elements of integer type T.
Iterator over points in a Mesh<D>.
void setDimensions(const IntVec< D > &dimensions)
Set the grid dimensions in all directions.
Exception thrown when not-a-number (NaN) is encountered.
Field of real double precision values on an FFT mesh.
System< D > const & system() const
Return const reference to parent system.
Compressor(System< D > &system)
Constructor.
int mdeCounter_
Count how many times MDE has been solved.
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.
LrCompressor(System< D > &system)
Constructor.
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.
void clearTimers()
Clear timers.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setup()
Initialize just before entry to iterative loop.
void outputTimers(std::ostream &out) const
Return compressor times contributions.
Main class, representing one complete 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.
Real periodic fields, SCFT and PS-FTS (CPU).
PSCF package top-level namespace.
float product(float a, float b)
Product for float Data.