1#ifndef RPG_LR_AM_COMPRESSOR_TPP
2#define RPG_LR_AM_COMPRESSOR_TPP
11#include "LrAmCompressor.h"
12#include <rpg/system/System.h>
13#include <rpg/solvers/Mixture.h>
14#include <rpg/field/Domain.h>
15#include <prdc/cuda/FFT.h>
16#include <prdc/cuda/resources.h>
17#include <pscf/mesh/MeshIterator.h>
32 isIntraCalculated_(false),
65 const int nMonomer = system().mixture().nMonomer();
66 const int meshSize = system().domain().mesh().size();
67 IntVec<D> const & dimensions = system().domain().mesh().dimensions();
75 for (
int i = 0; i < D; ++i) {
77 kMeshDimensions_[i] = dimensions[i];
79 kMeshDimensions_[i] = dimensions[i]/2 + 1;
85 for (
int i = 0; i < D; ++i) {
86 kSize_ *= kMeshDimensions_[i];
91 w0_.allocate(nMonomer);
92 wFieldTmp_.allocate(nMonomer);
93 resid_.allocate(dimensions);
94 residK_.allocate(dimensions);
95 intraCorrelationK_.allocate(kMeshDimensions_);
96 for (
int i = 0; i < nMonomer; ++i) {
97 w0_[i].allocate(dimensions);
98 wFieldTmp_[i].allocate(dimensions);
105 if (!isIntraCalculated_){
106 intra_.computeIntraCorrelations(intraCorrelationK_);
107 isIntraCalculated_ =
true;
111 for (
int i = 0; i < nMonomer; ++i) {
131 LrAmCompressor<D>::addCorrection(
133 VectorT
const & resTrial)
135 int n = fieldTrial.capacity();
136 const double vMonomer = system().mixture().vMonomer();
142 system().domain().fft().forwardTransform(resid_, residK_);
148 system().domain().fft().inverseTransformUnsafe(residK_, resid_);
160 int LrAmCompressor<D>::nElements()
161 {
return system().domain().mesh().size(); }
167 bool LrAmCompressor<D>::hasInitialGuess()
168 {
return system().w().hasData(); }
174 void LrAmCompressor<D>::getCurrent(VectorT& curr)
189 void LrAmCompressor<D>::evaluate()
199 void LrAmCompressor<D>::getResidual(VectorT& resid)
201 const int nMonomer = system().mixture().nMonomer();
207 for (
int i = 1; i < nMonomer; i++) {
216 void LrAmCompressor<D>::update(VectorT& newGuess)
219 const int nMonomer = system().mixture().nMonomer();
222 for (
int i = 0; i < nMonomer; i++) {
227 system().w().setRGrid(wFieldTmp_);
234 void LrAmCompressor<D>::outputToLog()
244 out <<
"LrAmCompressor time contributions:\n";
264 void LrAmCompressor<D>::setEqual(VectorT& a,
275 double LrAmCompressor<D>::dotProduct(VectorT
const & a,
286 double LrAmCompressor<D>::maxAbs(VectorT
const & a)
293 void LrAmCompressor<D>::subVV(VectorT& a,
302 void LrAmCompressor<D>::addEqVc(VectorT& a,
virtual void setup(bool isContinuation)
void outputTimers(std::ostream &out) const
void readErrorType(std::istream &in)
virtual void readParameters(std::istream &in)
int solve(bool isContinuation=false)
An IntVec<D, T> is a D-component vector of elements of integer type T.
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.
Base class for iterators that impose incompressibility.
int mdeCounter_
Count how many times MDE has been solved.
~LrAmCompressor()
Destructor.
void outputTimers(std::ostream &out) const
Return compressor times contributions.
void readParameters(std::istream &in)
Read all parameters and initialize.
int compress()
Compress to obtain partial saddle point w+.
void clearTimers()
Clear all timers (reset accumulated time to zero).
void setup(bool isContinuation)
Initialize just before entry to iterative loop.
LrAmCompressor(System< D > &system)
Constructor.
Main class, representing a complete physical system.
void setClassName(const char *className)
Set class name string.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
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 subVV(Array< double > &a, Array< double > const &b, Array< double > const &c)
Vector subtraction, a[i] = b[i] - c[i].
void addEqV(Array< double > &a, Array< double > const &b)
Vector addition in-place, a[i] += b[i].
void eqV(Array< double > &a, Array< double > const &b)
Vector assignment, 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.
void addEqVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c)
Vector addition in-place w/ coefficient, a[i] += b[i] * c, kernel wrapper.
SCFT and PS-FTS with real periodic fields (GPU)
PSCF package top-level namespace.