1#ifndef RPC_LR_POST_AM_COMPRESSOR_TPP
2#define RPC_LR_POST_AM_COMPRESSOR_TPP
11#include "LrAmCompressor.h"
12#include <rpc/System.h>
13#include <rpc/fts/compressor/intra/IntraCorrelation.h>
14#include <pscf/mesh/MeshIterator.h>
27 isIntraCalculated_(false),
29 { setClassName(
"LrAmCompressor"); }
54 const int nMonomer = system().mixture().nMonomer();
55 const int meshSize = system().domain().mesh().size();
56 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
59 for (
int i = 0; i < D; ++i) {
61 kMeshDimensions_[i] = dimensions[i];
63 kMeshDimensions_[i] = dimensions[i]/2 + 1;
69 newBasis_.allocate(meshSize);
70 w0_.allocate(nMonomer);
71 wFieldTmp_.allocate(nMonomer);
72 resid_.allocate(dimensions);
73 residK_.allocate(dimensions);
74 intraCorrelationK_.allocate(kMeshDimensions_);
75 for (
int i = 0; i < nMonomer; ++i) {
76 w0_[i].allocate(dimensions);
77 wFieldTmp_[i].allocate(dimensions);
84 if (!isIntraCalculated_){
85 intra_.computeIntraCorrelations(intraCorrelationK_);
86 isIntraCalculated_ =
true;
90 for (
int i = 0; i < nMonomer; ++i) {
91 for (
int j = 0; j< meshSize; ++j){
92 w0_[i][j] = system().w().rgrid(i)[j];
127 for (
int i = 0; i < n; i++) {
129 if (std::isnan(a[i]) || std::isnan(b[i])) {
131 __FILE__,__LINE__,0);
147 for (
int i = 0; i < n; i++) {
149 if (std::isnan(value)) {
150 throw NanException(
"LrAmCompressor::dotProduct",
151 __FILE__, __LINE__, 0);
153 if (fabs(value) > max) {
165 LrAmCompressor<D>::updateBasis(
172 const int n = hists[0].capacity();
175 for (
int i = 0; i < n; i++) {
176 newBasis_[i] = hists[0][i] - hists[1][i];
179 basis.append(newBasis_);
190 for (
int i = 0; i < nHist; i++) {
191 for (
int j = 0; j < n; j++) {
193 trial[j] += coeffs[i] * -1 * basis[i][j];
202 double LrAmCompressor<D>::computeLambda(
double r)
218 const int n = fieldTrial.
capacity();
219 const double vMonomer = system().mixture().vMonomer();
223 for (
int i = 0 ; i < n; ++i) {
224 resid_[i] = resTrial[i];
228 system().domain().fft().forwardTransform(resid_, residK_);
232 MeshIterator<D> iter;
233 iter.setDimensions(kMeshDimensions_);
234 for (iter.begin(); !iter.atEnd(); ++iter) {
235 factor = 1.0 / (vMonomer * intraCorrelationK_[iter.rank()]);
236 residK_[iter.rank()][0] *= factor;
237 residK_[iter.rank()][1] *= factor;
243 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
246 for (
int i = 0; i < n; i++) {
247 fieldTrial[i] += lambda * resid_[i];
255 bool LrAmCompressor<D>::hasInitialGuess()
256 {
return system().w().hasData(); }
262 int LrAmCompressor<D>::nElements()
263 {
return system().domain().mesh().size(); }
275 const int meshSize = system().domain().mesh().size();
284 for (
int i = 0; i < meshSize; i++){
285 curr[i] = (*currSys)[0][i] - w0_[0][i];
294 void LrAmCompressor<D>::evaluate()
306 const int n = nElements();
307 const int nMonomer = system().mixture().nMonomer();
308 const int meshSize = system().domain().mesh().size();
311 for (
int i = 0 ; i < n; ++i) {
316 for (
int j = 0; j < nMonomer; ++j) {
317 for (
int k = 0; k < meshSize; ++k) {
318 resid[k] += system().c().rgrid(j)[k];
331 const int nMonomer = system().mixture().nMonomer();
332 const int meshSize = system().domain().mesh().size();
336 for (
int i = 0; i < nMonomer; i++){
337 for (
int k = 0; k < meshSize; k++){
338 wFieldTmp_[i][k] = w0_[i][k] + newGuess[k];
343 system().setWRGrid(wFieldTmp_);
347 void LrAmCompressor<D>::outputToLog()
358 out <<
"LrAmCompressor time contributions:\n";
Template for Anderson mixing iterator algorithm.
An IntVec<D, T> is a D-component vector of elements of integer type T.
Exception thrown when not-a-number (NaN) is encountered.
Base class for iterators that impose incompressibility.
Anderson Mixing compressor with linear-response mixing step.
void readParameters(std::istream &in)
Read all parameters and initialize.
~LrAmCompressor()
Destructor.
LrAmCompressor(System< D > &system)
Constructor.
void outputTimers(std::ostream &out)
Return compressor times contributions.
void clearTimers()
Clear all timers (reset accumulated time to zero).
void setup(bool isContinuation)
Initialize just before entry to iterative loop.
int compress()
Compress to obtain partial saddle point w+.
Main class for SCFT or PS-FTS simulation of one system.
int capacity() const
Return allocated size.
Dynamically allocatable contiguous array template.
Class for storing history of previous values in an array.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
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.