1#ifndef RPG_AM_ITERATOR_GRID_TPP
2#define RPG_AM_ITERATOR_GRID_TPP
11#include "AmIteratorGrid.h"
12#include <rpg/System.h>
13#include <prdc/crystal/UnitCell.h>
14#include <prdc/cuda/resources.h>
15#include <pscf/inter/Interaction.h>
16#include <pscf/iterator/NanException.h>
26 using namespace Prdc::Cuda;
32 imposedFields_(system)
52 interaction_.setNMonomer(system().mixture().nMonomer());
58 int np = system().domain().unitCell().nParameter();
64 readOptional(in,
"isFlexible", isFlexible_);
69 flexibleParams_.clear();
70 for (
int i = 0; i < np; i++) {
71 flexibleParams_.append(
true);
74 readOptionalFSArray(in,
"flexibleParams", flexibleParams_, np);
75 if (nFlexibleParams() == 0) isFlexible_ =
false;
77 flexibleParams_.clear();
78 for (
int i = 0; i < np; i++) {
79 flexibleParams_.append(
false);
84 readOptional(in,
"scaleStress", scaleStress_);
87 readParamCompositeOptional(in, imposedFields_);
96 out <<
"Iterator times contributions:\n";
106 if (imposedFields_.isActive()) {
107 imposedFields_.setup();
111 interaction_.update(system().interaction());
126 double AmIteratorGrid<D>::dotProduct(FieldCUDA
const & a,
131 if (std::isnan(val)) {
132 throw NanException(
"AmIteratorGrid::dotProduct", __FILE__,
140 double AmIteratorGrid<D>::maxAbs(FieldCUDA
const & a)
143 if (std::isnan(val)) {
144 throw NanException(
"AmIteratorGrid::maxAbs", __FILE__, __LINE__, 0);
159 if (basis[0].isAllocated()) {
160 UTIL_CHECK(basis[0].capacity() == hists[0].capacity());
162 basis[0].
allocate(hists[0].capacity());
172 AmIteratorGrid<D>::addHistories(FieldCUDA& trial,
176 for (
int i = 0; i < nHist; i++) {
183 void AmIteratorGrid<D>::addPredictedError(FieldCUDA& fieldTrial,
184 FieldCUDA
const & resTrial,
192 bool AmIteratorGrid<D>::hasInitialGuess()
193 {
return system().w().hasData(); }
197 int AmIteratorGrid<D>::nElements()
199 const int nMonomer = system().mixture().nMonomer();
200 const int nMesh = system().mesh().size();
202 int nEle = nMonomer*nMesh;
205 nEle += nFlexibleParams();
213 void AmIteratorGrid<D>::getCurrent(FieldCUDA& curr)
215 const int nMonomer = system().mixture().nMonomer();
216 const int nMesh = system().mesh().size();
217 const int n = nElements();
221 for (
int i = 0; i < nMonomer; i++) {
222 VecOp::eqV(curr, system().w().rgrid(i), i*nMesh, 0, nMesh);
227 const int nParam = system().unitCell().nParameter();
229 = system().unitCell().parameters();
232 HostDArray<cudaReal> tempH(nFlexibleParams());
234 for (
int i = 0; i < nParam; i++) {
235 if (flexibleParams_[i]) {
236 tempH[counter] = scaleStress_ * parameters[i];
244 tempD.associate(curr, nMonomer*nMesh, tempH.capacity());
251 void AmIteratorGrid<D>::evaluate()
254 system().compute(isFlexible_);
259 void AmIteratorGrid<D>::getResidual(FieldCUDA& resid)
261 const int n = nElements();
262 const int nMonomer = system().mixture().nMonomer();
263 const int nMesh = system().mesh().size();
273 for (
int i = 0; i < nMonomer; i++) {
274 residSlices[i].associate(resid, i*nMesh, nMesh);
278 for (
int i = 0; i < nMonomer; i++) {
279 for (
int j = 0; j < nMonomer; j++) {
281 system().c().rgrid(j), interaction_.chi(i, j),
282 system().w().rgrid(j), -interaction_.p(i, j));
287 if (system().hasMask()) {
288 double coeff = -1.0 / interaction_.sumChiInverse();
289 for (
int i = 0; i < nMonomer; ++i) {
296 if (system().hasExternalFields()) {
297 for (
int i = 0; i < nMonomer; ++i) {
298 for (
int j = 0; j < nMonomer; ++j) {
299 double p = interaction_.p(i,j);
306 if (!system().mixture().isCanonical()) {
307 cudaReal factor = 1.0 / interaction_.sumChiInverse();
310 for (
int i = 0; i < nMonomer; i++) {
312 cudaReal average = findAverage(residSlices[i]);
320 const int nParam = system().unitCell().nParameter();
321 HostDArray<cudaReal> stressH(nFlexibleParams());
324 for (
int i = 0; i < nParam; i++) {
325 if (flexibleParams_[i]) {
326 double stress = system().mixture().stress(i);
329 if (imposedFields_.isActive()) {
330 stress = imposedFields_.correctedStress(i,stress);
333 stressH[counter] = -1 * scaleStress_ * stress;
340 stressD.associate(resid, nMonomer*nMesh, stressH.capacity());
347 void AmIteratorGrid<D>::update(FieldCUDA& newGuess)
349 const int nMonomer = system().mixture().nMonomer();
350 const int nMesh = system().mesh().size();
353 if (system().mixture().isCanonical()) {
354 cudaReal average, wAverage, cAverage;
355 for (
int i = 0; i < nMonomer; i++) {
358 ngSlice.associate(newGuess, i*nMesh, nMesh);
361 average = findAverage(ngSlice);
368 for (
int j = 0; j < nMonomer; j++) {
370 if (system().w().isSymmetric()) {
371 cAverage = system().c().basis(j)[0];
373 cAverage = findAverage(system().c().rgrid(j));
375 wAverage += interaction_.chi(i,j) * cAverage;
379 if (system().hasExternalFields()) {
380 if (system().h().isSymmetric()) {
381 wAverage += system().h().basis(i)[0];
383 wAverage += findAverage(system().h().rgrid(i));
392 system().setWRGrid(newGuess);
393 system().symmetrizeWFields();
398 parameters = system().domain().unitCell().parameters();
399 const int nParam = system().unitCell().nParameter();
402 HostDArray<cudaReal> tempH(nFlexibleParams());
403 tempH.copySlice(newGuess, nMonomer*nMesh);
405 const double coeff = 1.0 / scaleStress_;
407 for (
int i = 0; i < nParam; i++) {
408 if (flexibleParams_[i]) {
409 parameters[i] = coeff * tempH[i];
415 system().setUnitCell(parameters);
419 if (imposedFields_.isActive()) {
420 imposedFields_.update();
426 void AmIteratorGrid<D>::outputToLog()
428 if (isFlexible_ && verbose() > 1) {
429 const int nParam = system().domain().unitCell().nParameter();
430 const int nMonomer = system().mixture().nMonomer();
431 const int nMesh = system().mesh().size();
434 HostDArray<cudaReal> tempH(nFlexibleParams());
435 tempH.copySlice(residual(), nMonomer*nMesh);
438 for (
int i = 0; i < nParam; i++) {
439 if (flexibleParams_[i]) {
440 double stress = tempH[counter] / (-1.0 * scaleStress_);
442 <<
" Cell Param " << i <<
" = "
443 <<
Dbl(system().domain().unitCell().parameters()[i], 15)
458 if (imposedFields_.isActive()) {
459 arr = imposedFields_.getParameterTypes();
466 double value,
bool& success)
468 if (imposedFields_.isActive()) {
469 imposedFields_.setParameter(name, ids, value, success);
480 if (imposedFields_.isActive()) {
481 return imposedFields_.getParameter(name, ids, success);
Template for Anderson mixing iterator algorithm.
Dynamic array on the GPU device with aligned data.
int capacity() const
Return allocated capacity.
Exception thrown when not-a-number (NaN) is encountered.
Base template for UnitCell<D> classes, D=1, 2 or 3.
Rpg implementation of the Anderson Mixing iterator.
AmIteratorGrid(System< D > &system)
Constructor.
void readParameters(std::istream &in)
Read all parameters and initialize.
void setClassName(const char *className)
Set class name string.
void setParameter(std::string name, DArray< int > ids, double value, bool &success)
Set the value of a specialized sweep parameter.
double getParameter(std::string name, DArray< int > ids, bool &success) const
Get the value of a specialized sweep parameter.
void setup(bool isContinuation)
Setup iterator just before entering iteration loop.
void outputTimers(std::ostream &out)
Output timing results to log file.
GArray< ParameterType > getParameterTypes()
Return specialized sweep parameter types to add to the Sweep object.
~AmIteratorGrid()
Destructor.
Base class for iterative solvers for SCF equations.
Main class for calculations that represent one system.
Dynamically allocatable contiguous array template.
void allocate(int capacity)
Allocate the underlying C array.
Wrapper for a double precision number, for formatted ostream output.
A fixed capacity (static) contiguous array with a variable logical size.
An automatically growable array, analogous to a std::vector.
static std::ostream & file()
Get log ostream by reference.
Class for storing history of previous values in an array.
void allocate(int capacity)
Allocate a new empty buffer.
void advance()
Advances the pointer to an added element without assigning a value.
int size() const
Return number of values currently in the buffer.
File containing preprocessor macros for error handling.
#define UTIL_CHECK(condition)
Assertion macro suitable for serial or parallel production code.
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).
cudaReal sum(DeviceArray< cudaReal > const &in)
Compute sum of array elements (GPU kernel wrapper).
void eqV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, const int beginIdA, const int beginIdB, const int n)
Vector assignment, a[i] = b[i], kernel wrapper (cudaReal).
void subEqS(DeviceArray< cudaReal > &a, cudaReal const b, const int beginIdA, const int n)
Vector subtraction in-place, a[i] -= b, kernel wrapper (cudaReal).
void addEqS(DeviceArray< cudaReal > &a, cudaReal const b, const int beginIdA, const int n)
Vector addition in-place, a[i] += b, kernel wrapper (cudaReal).
void addVcVcVc(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, cudaReal const c, DeviceArray< cudaReal > const &d, cudaReal const e, DeviceArray< cudaReal > const &f, cudaReal const g)
3-vec addition w coeff, a[i] = (b[i]*c) + (d[i]*e) + (f[i]*g), kernel wrapper.
void subVV(DeviceArray< cudaReal > &a, DeviceArray< cudaReal > const &b, DeviceArray< cudaReal > const &c, const int beginIdA, const int beginIdB, const int beginIdC, const int n)
Vector subtraction, a[i] = b[i] - c[i], kernel wrapper (cudaReal).
void eqS(DeviceArray< cudaReal > &a, cudaReal const b, const int beginIdA, const int n)
Vector assignment, a[i] = b, kernel wrapper (cudaReal).
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.
PSCF package top-level namespace.
Utility classes for scientific computation.