← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/atc/Kinetostat.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2013-11-20 22:41:36 UTC
  • mfrom: (1.2.2)
  • Revision ID: package-import@ubuntu.com-20131120224136-tzx7leh606fqnckm
Tags: 0~20131119.git7162cf0-1
* [e65b919] Imported Upstream version 0~20131119.git7162cf0
* [f7bddd4] Fix some problems, introduced by upstream recently.
* [3616dfc] Use wrap-and-sort script.
* [7e92030] Ignore quilt dir

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/atc/Kinetostat.cpp
1
 
// ATC_Transfer Headers
 
1
#include "Kinetostat.h"
2
2
#include "ATC_Error.h"
3
 
#include "Kinetostat.h"
4
 
#include "CG.h"
5
 
#include "ATC_Transfer.h"
 
3
#include "ATC_Coupling.h"
6
4
#include "LammpsInterface.h"
 
5
#include "PerAtomQuantityLibrary.h"
 
6
#include "PrescribedDataManager.h"
 
7
#include "ElasticTimeIntegrator.h"
 
8
#include "TransferOperator.h"
 
9
 
 
10
using std::set;
 
11
using std::pair;
 
12
using std::string;
7
13
 
8
14
namespace ATC {
9
15
 
16
22
  //--------------------------------------------------------
17
23
  //  Constructor
18
24
  //--------------------------------------------------------
19
 
  Kinetostat::Kinetostat(ATC_Transfer * atcTransfer) :
20
 
    AtomicRegulator(atcTransfer),
21
 
    kinetostatType_(NONE),
22
 
    couplingMode_(UNCOUPLED)
 
25
  Kinetostat::Kinetostat(ATC_Coupling * atc,
 
26
                         const string & regulatorPrefix) :
 
27
    AtomicRegulator(atc,regulatorPrefix)
23
28
  {
24
29
    // do nothing
25
30
  }
33
38
  {
34
39
    bool foundMatch = false;
35
40
 
36
 
    // fluxstat type
37
 
    /*! \page man_disp_control fix_modify AtC transfer momentum control
38
 
      \section syntax
39
 
      fix_modify AtC transfer momentum control none \n
40
 
      
41
 
      fix_modify AtC transfer momentum control rescale <frequency>\n
42
 
      - frequency (int) = time step frequency for applying displacement and velocity rescaling \n
43
 
      
44
 
      fix_modify AtC transfer momentum control glc_displacement \n
45
 
  
46
 
      fix_modify AtC transfer momentum control glc_velocity \n
47
 
 
48
 
      fix_modify AtC transfer momentum control flux [faceset face_set_id, interpolate] 
49
 
      - face_set_id (string) = id of boundary face set, if not specified
50
 
      (or not possible when the atomic domain does not line up with 
51
 
      mesh boundaries) defaults to an atomic-quadrature approximate 
52
 
      evaulation\n
53
 
      \section examples
54
 
      fix_modify AtC transfer momentum control glc_velocity \n
55
 
      fix_modify AtC transfer momentum control stress_flux faceset bndy_faces \n
56
 
      \section description
57
 
      \section restrictions
58
 
      only for be used with specific transfers :
59
 
      elastic \n
60
 
      rescale not valid with time filtering activated
61
 
      \section related
62
 
      \section default
63
 
      none
64
 
    */
65
 
    
66
41
    int argIndex = 0;
67
 
    if (strcmp(arg[argIndex],"none")==0) { // restore defaults
68
 
      kinetostatType_ = NONE;
69
 
      couplingMode_ = UNCOUPLED;
70
 
      howOften_ = 1;
71
 
      boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
72
 
      foundMatch = true;
73
 
    }
74
 
    else if (strcmp(arg[argIndex],"glc_displacement")==0) {
75
 
      kinetostatType_ = GLC_DISPLACEMENT;
76
 
      couplingMode_ = FIXED;
77
 
      howOften_ = 1;
78
 
      boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
79
 
      foundMatch = true;
80
 
    }
81
 
    else if (strcmp(arg[argIndex],"glc_velocity")==0) {
82
 
      kinetostatType_ = GLC_VELOCITY;
83
 
      couplingMode_ = FIXED;
84
 
      howOften_ = 1;
85
 
      boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
86
 
      foundMatch = true;
87
 
    }
88
 
    else if (strcmp(arg[argIndex],"hoover")==0) {
89
 
      kinetostatType_ = FORCE;
90
 
      couplingMode_ = FIXED;
91
 
      howOften_ = 1;
92
 
      boundaryIntegrationType_ = ATC_Transfer::NO_QUADRATURE;
93
 
      foundMatch = true;
94
 
    }
95
 
    else if (strcmp(arg[argIndex],"flux")==0) {
96
 
      kinetostatType_ = FORCE;
97
 
      couplingMode_ = FLUX;
98
 
      howOften_ = 1;
 
42
    if (strcmp(arg[argIndex],"momentum")==0) {
99
43
      argIndex++;
100
 
      boundaryIntegrationType_ = atcTransfer_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_);
101
 
      foundMatch = true;
 
44
 
 
45
      // fluxstat type
 
46
      /*! \page man_control_momentum fix_modify AtC control momentum
 
47
        \section syntax
 
48
        fix_modify AtC control momentum none \n
 
49
      
 
50
        fix_modify AtC control momentum rescale <frequency>\n
 
51
        - frequency (int) = time step frequency for applying displacement and velocity rescaling \n
 
52
      
 
53
        fix_modify AtC control momentum glc_displacement \n
 
54
  
 
55
        fix_modify AtC control momentum glc_velocity \n
 
56
 
 
57
        fix_modify AtC control momentum hoover \n
 
58
 
 
59
        fix_modify AtC control momentum flux [faceset face_set_id, interpolate] 
 
60
        - face_set_id (string) = id of boundary face set, if not specified
 
61
        (or not possible when the atomic domain does not line up with 
 
62
        mesh boundaries) defaults to an atomic-quadrature approximate 
 
63
        evaulation\n
 
64
        \section examples
 
65
        fix_modify AtC control momentum glc_velocity \n
 
66
        fix_modify AtC control momentum flux faceset bndy_faces \n
 
67
        \section description
 
68
        \section restrictions
 
69
        only to be used with specific transfers :
 
70
        elastic \n
 
71
        rescale not valid with time filtering activated
 
72
        \section related
 
73
        \section default
 
74
        none
 
75
      */
 
76
      boundaryIntegrationType_ = NO_QUADRATURE;
 
77
      howOften_ = 1;
 
78
      if (strcmp(arg[argIndex],"none")==0) { // restore defaults
 
79
        regulatorTarget_ = NONE;
 
80
        couplingMode_ = UNCOUPLED;
 
81
        foundMatch = true;
 
82
      }
 
83
      else if (strcmp(arg[argIndex],"glc_displacement")==0) {
 
84
        regulatorTarget_ = FIELD;
 
85
        couplingMode_ = FIXED;
 
86
        foundMatch = true;
 
87
      }
 
88
      else if (strcmp(arg[argIndex],"glc_velocity")==0) {
 
89
        regulatorTarget_ = DERIVATIVE;
 
90
        couplingMode_ = FIXED;
 
91
        foundMatch = true;
 
92
      }
 
93
      else if (strcmp(arg[argIndex],"hoover")==0) {
 
94
        regulatorTarget_ = DYNAMICS;
 
95
        couplingMode_ = FIXED;
 
96
        foundMatch = true;
 
97
      }
 
98
      else if (strcmp(arg[argIndex],"flux")==0) {
 
99
        regulatorTarget_ = DYNAMICS;
 
100
        couplingMode_ = FLUX;
 
101
        argIndex++;
 
102
        boundaryIntegrationType_ = atc_->parse_boundary_integration(narg-argIndex,&arg[argIndex],boundaryFaceSet_);
 
103
        foundMatch = true;
 
104
      }
 
105
      else if (strcmp(arg[argIndex],"ghost_flux")==0) {
 
106
        regulatorTarget_ = DYNAMICS;
 
107
        couplingMode_ = GHOST_FLUX;
 
108
        foundMatch = true;
 
109
      }
102
110
    }
103
111
  
104
112
    if (!foundMatch)
109
117
  }
110
118
 
111
119
  //--------------------------------------------------------
112
 
  //  reset_nodal_lambda_force:
 
120
  //  reset_lambda_contribution
113
121
  //    resets the kinetostat generated force to a
114
122
  //    prescribed value
115
123
  //--------------------------------------------------------
116
 
  void Kinetostat::reset_lambda_force(DENS_MAT & target)
 
124
  void Kinetostat::reset_lambda_contribution(const DENS_MAT & target)
117
125
  {
118
 
    for (int i = 0; i < nNodes_; ++i)
119
 
      for (int j = 0; j < nsd_; ++j)
120
 
        lambdaForceFiltered_(i,j) = target(i,j);
 
126
    DENS_MAN * lambdaForceFiltered = regulator_data("LambdaForceFiltered",nsd_);
 
127
    lambdaForceFiltered->set_quantity() = target;
121
128
  }
122
129
 
123
130
  //--------------------------------------------------------
124
131
  //  initialize:
125
132
  //    sets up methods before a run
126
 
  //    NOTE we currently only include time filter
 
133
  
127
134
  //    dependence, but in general there is also a
128
135
  //    time integrator dependence.  In general the 
129
136
  //    precedence order is:
132
139
  //    different types of time integrators can be
133
140
  //    specified.
134
141
  //--------------------------------------------------------
135
 
  void Kinetostat::initialize()
 
142
  void Kinetostat::construct_methods()
136
143
  {
137
 
    // NOTE: only support basic time integration
138
 
    TimeFilterManager * timeFilterManager = atcTransfer_->get_time_filter_manager();
139
 
 
140
 
    // reset data if needed and perform any error/conflict checking
141
 
    if (resetData_) {
142
 
      AtomicRegulator::reset_data();
143
 
 
144
 
      // set up storage
145
 
      lambda_.reset(nNodes_,nsd_);
146
 
      nodalAtomicLambdaForce_.reset(nNodes_,nsd_);
147
 
      lambdaForceFiltered_.reset(nNodes_,nsd_);
148
 
    }
149
 
    
150
 
    if (needReset_ || timeFilterManager->need_reset() || timeFilterManager->end_equilibrate()) {
 
144
    // get data associated with stages 1 & 2 of ATC_Method::initialize
 
145
    AtomicRegulator::construct_methods();
 
146
 
 
147
    if (atc_->reset_methods()) {
151
148
      // eliminate existing methods
152
 
      destroy();
 
149
      delete_method();
153
150
      
154
151
      DENS_MAT nodalGhostForceFiltered;
155
 
      if (timeFilterManager->end_equilibrate() && kinetostatType_==FORCE) {
 
152
      TimeIntegrator::TimeIntegrationType myIntegrationType = (atc_->time_integrator(VELOCITY))->time_integration_type();
 
153
      TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
 
154
      if (timeFilterManager->end_equilibrate() && regulatorTarget_==AtomicRegulator::DYNAMICS) {
156
155
        StressFlux * myMethod;
157
156
        myMethod = dynamic_cast<StressFlux *>(regulatorMethod_);
158
 
        nodalGhostForceFiltered = myMethod->get_filtered_ghost_force();
 
157
        nodalGhostForceFiltered = (myMethod->filtered_ghost_force()).quantity();
159
158
      }
160
 
 
 
159
      
 
160
      // update time filter
 
161
      
161
162
      if (timeFilterManager->need_reset()) {
162
 
        if (timeFilter_)
163
 
          delete timeFilter_;
 
163
        
164
164
        timeFilter_ = timeFilterManager->construct(TimeFilterManager::IMPLICIT_UPDATE);
165
165
      }
166
166
                         
167
167
      if (timeFilterManager->filter_dynamics()) {
168
 
        switch (kinetostatType_) {
 
168
        switch (regulatorTarget_) {
169
169
        case NONE: {
 
170
          regulatorMethod_ = new RegulatorMethod(this);
170
171
          break;
171
172
        }
172
 
        case GLC_DISPLACEMENT: {
 
173
        case FIELD: {
173
174
          regulatorMethod_ = new DisplacementGlcFiltered(this);
174
175
          break;
175
176
        }
176
 
        case GLC_VELOCITY: {
 
177
        case DERIVATIVE: {
177
178
          regulatorMethod_ = new VelocityGlcFiltered(this);
178
179
          break;
179
180
        }
180
 
        case FORCE: {
 
181
        case DYNAMICS: {
 
182
          throw ATC_Error("Kinetostat::initialize - force based kinetostats not yet implemented with time filtering");
181
183
          regulatorMethod_ = new StressFluxFiltered(this);
182
184
          if (timeFilterManager->end_equilibrate()) {
183
185
            StressFlux * myMethod;
187
189
          break;
188
190
        }
189
191
        default:
190
 
          throw ATC_Error(0,"Unknown kinetostat type in Kinetostat::initialize");
 
192
          throw ATC_Error("Unknown kinetostat type in Kinetostat::initialize");
191
193
        }
192
194
      }
193
195
      else {
194
 
        switch (kinetostatType_) {
 
196
        switch (regulatorTarget_) {
195
197
        case NONE: {
 
198
          regulatorMethod_ = new RegulatorMethod(this);
196
199
          break;
197
200
        }
198
 
        case GLC_DISPLACEMENT: {
 
201
        case FIELD: {
199
202
          regulatorMethod_ = new DisplacementGlc(this);
200
203
          break;
201
204
        }
202
 
        case GLC_VELOCITY: {
 
205
        case DERIVATIVE: {
203
206
          regulatorMethod_ = new VelocityGlc(this);
204
207
          break;
205
208
        }
206
 
        case FORCE: {
207
 
          regulatorMethod_ = new StressFlux(this);
208
 
          break;
 
209
        case DYNAMICS: {
 
210
          if (myIntegrationType == TimeIntegrator::FRACTIONAL_STEP) {
 
211
            if (couplingMode_ == GHOST_FLUX) {
 
212
              regulatorMethod_ = new KinetostatFluxGhost(this);
 
213
            }
 
214
            else if (couplingMode_ == FIXED) {
 
215
              regulatorMethod_ = new KinetostatFixed(this);
 
216
            }
 
217
            else if (couplingMode_ == FLUX) {
 
218
              regulatorMethod_ = new KinetostatFlux(this);
 
219
            }
 
220
            break;
 
221
          }
 
222
          if (myIntegrationType == TimeIntegrator::GEAR) {
 
223
            if (couplingMode_ == FIXED) {
 
224
              regulatorMethod_ = new KinetostatFixed(this);
 
225
            }
 
226
            else if (couplingMode_ == FLUX) {
 
227
              regulatorMethod_ = new KinetostatFlux(this);
 
228
            }
 
229
            break;
 
230
          }
 
231
          else {
 
232
            if (couplingMode_ == GHOST_FLUX) {
 
233
              regulatorMethod_ = new StressFluxGhost(this);
 
234
            }
 
235
            else {
 
236
              regulatorMethod_ = new StressFlux(this);
 
237
            }
 
238
            break;
 
239
          }
209
240
        }
210
241
        default:
211
 
          throw ATC_Error(0,"Unknown kinetostat type in Kinetostat::initialize");
 
242
          throw ATC_Error("Unknown kinetostat type in Kinetostat::initialize");
212
243
        }
 
244
        AtomicRegulator::reset_method();
213
245
      }
214
 
 
215
 
      AtomicRegulator::reset_method();
216
 
    }
217
 
 
218
 
     AtomicRegulator::initialize();
 
246
    }
 
247
    else {
 
248
      set_all_data_to_used();
 
249
    }
219
250
  }
 
251
 
 
252
  //--------------------------------------------------------
 
253
  //--------------------------------------------------------
 
254
  //  Class KinetostatShapeFunction
 
255
  //--------------------------------------------------------
 
256
  //--------------------------------------------------------
220
257
 
221
258
  //--------------------------------------------------------
 
259
  //  Constructor
 
260
  //         Grab references to ATC and kinetostat data
 
261
  //--------------------------------------------------------
 
262
  KinetostatShapeFunction::KinetostatShapeFunction(Kinetostat *kinetostat,
 
263
                                                   const string & regulatorPrefix) :
 
264
    RegulatorShapeFunction(kinetostat,regulatorPrefix),
 
265
    kinetostat_(kinetostat),
 
266
    timeFilter_(atomicRegulator_->time_filter()),
 
267
    nodalAtomicLambdaForce_(NULL),
 
268
    lambdaForceFiltered_(NULL),
 
269
    atomKinetostatForce_(NULL),
 
270
    atomVelocities_(NULL),
 
271
    atomMasses_(NULL)
 
272
  {
 
273
    // data associated with stage 3 in ATC_Method::initialize
 
274
    lambda_ = kinetostat->regulator_data(regulatorPrefix_+"LambdaMomentum",nsd_);
 
275
    lambdaForceFiltered_ = kinetostat_->regulator_data("LambdaForceFiltered",nsd_);
 
276
  }
 
277
 
 
278
  //--------------------------------------------------------
 
279
  //  constructor_transfers
 
280
  //    instantiates or obtains all dependency managed data
 
281
  //--------------------------------------------------------
 
282
  void KinetostatShapeFunction::construct_transfers()
 
283
  {
 
284
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
285
 
 
286
    // needed fundamental quantities
 
287
    atomVelocities_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_VELOCITY);
 
288
    atomMasses_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_MASS);
 
289
 
 
290
    // base class transfers
 
291
    RegulatorShapeFunction::construct_transfers();
 
292
 
 
293
    // lambda interpolated to the atomic coordinates
 
294
    atomLambdas_ = new FtaShapeFunctionProlongation(atc_,
 
295
                                                   lambda_,
 
296
                                                    interscaleManager.per_atom_sparse_matrix("Interpolant"));
 
297
    interscaleManager.add_per_atom_quantity(atomLambdas_,
 
298
                                            regulatorPrefix_+"AtomLambdaMomentum");
 
299
  }
 
300
 
 
301
  //--------------------------------------------------------
 
302
  //  set_weights
 
303
  //            sets diagonal weighting matrix used in
 
304
  //            solve_for_lambda
 
305
  //--------------------------------------------------------
 
306
  void KinetostatShapeFunction::set_weights()
 
307
  {
 
308
    if (this->use_local_shape_functions()) {
 
309
      ConstantQuantityMapped<double> * myWeights = new ConstantQuantityMapped<double>(atc_,1.,lambdaAtomMap_);
 
310
      weights_ = myWeights;
 
311
      (atc_->interscale_manager()).add_per_atom_quantity(myWeights,
 
312
                                                         "AtomOnesMapped");
 
313
    }
 
314
    else {
 
315
      weights_ = (atc_->interscale_manager()).per_atom_quantity("AtomicOnes");
 
316
      if (!weights_) {
 
317
        weights_ = new ConstantQuantity<double>(atc_,1.);
 
318
        (atc_->interscale_manager()).add_per_atom_quantity(weights_,
 
319
                                                           "AtomicOnes");
 
320
      }
 
321
    }
 
322
  }
 
323
 
 
324
  //--------------------------------------------------------
222
325
  //--------------------------------------------------------
223
326
  //  Class GlcKinetostat
224
327
  //--------------------------------------------------------
228
331
  //  Constructor
229
332
  //         Grab references to ATC and kinetostat data
230
333
  //--------------------------------------------------------
231
 
  GlcKinetostat::GlcKinetostat(Kinetostat *kinetostat) 
232
 
    :
233
 
    RegulatorShapeFunction(kinetostat),
234
 
    kinetostat_(kinetostat),
235
 
    timeFilter_(atomicRegulator_->get_time_filter()),
236
 
    nodalAtomicLambdaForce_(kinetostat_->get_nodal_atomic_lambda_force()),
237
 
    lambdaForceFiltered_(kinetostat_->get_lambda_force_filtered()),
238
 
    lambdaForce_(atomicRegulator_->get_lambda_force()),
239
 
    mdMassMatrix_(atcTransfer_->get_mass_mat_md(VELOCITY)),
240
 
    nodeToOverlapMap_(atcTransfer_->get_node_to_overlap_map())
241
 
  {
 
334
  GlcKinetostat::GlcKinetostat(Kinetostat *kinetostat) :
 
335
    KinetostatShapeFunction(kinetostat),
 
336
    mdMassMatrix_(atc_->set_mass_mat_md(VELOCITY)),
 
337
    atomPositions_(NULL)
 
338
  {
 
339
    // do nothing
 
340
  }
 
341
 
 
342
  //--------------------------------------------------------
 
343
  //  constructor_transfers
 
344
  //    instantiates or obtains all dependency managed data
 
345
  //--------------------------------------------------------
 
346
  void GlcKinetostat::construct_transfers()
 
347
  {
 
348
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
349
 
 
350
    // needed fundamental quantities
 
351
    atomPositions_ = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_POSITION);
 
352
    
 
353
    // base class transfers
 
354
    KinetostatShapeFunction::construct_transfers();
 
355
  }
 
356
 
 
357
  //--------------------------------------------------------
 
358
  //  initialize
 
359
  //    initializes all method data
 
360
  //--------------------------------------------------------
 
361
  void GlcKinetostat::initialize()
 
362
  {
 
363
    KinetostatShapeFunction::initialize();
 
364
 
 
365
    
242
366
    // set up list of nodes using Hoover coupling
243
367
    // (a) nodes with prescribed values
244
 
    Array2D<bool> isFixedNode = atcTransfer_->get_fixed_node_flags(VELOCITY);
 
368
    PrescribedDataManager * prescribedDataMgr(atc_->prescribed_data_manager());
245
369
    for (int i = 0; i < nNodes_; ++i)
246
 
      for (int j = 0; j < isFixedNode.nCols(); ++j)
247
 
        if (isFixedNode(i,j))
 
370
      for (int j = 0; j < nsd_; ++j)
 
371
        if (prescribedDataMgr->is_fixed(i,VELOCITY,j))
248
372
          hooverNodes_.insert(pair<int,int>(i,j));
249
373
 
250
374
    // (b) AtC coupling nodes
251
 
    if (kinetostat_->get_coupling_mode()==Kinetostat::FIXED) {
252
 
      Array<int> & nodeType(atcTransfer_->get_node_type());
253
 
      if (atcTransfer_->use_localized_lambda()) {
254
 
        for (int i = 0; i < nNodes_; ++i)
255
 
          if (nodeType(i)==ATC_Transfer::BOUNDARY)
256
 
            for (int j = 0; j < nsd_; ++j)
 
375
    if (atomicRegulator_->coupling_mode()==AtomicRegulator::FIXED) {
 
376
      InterscaleManager & interscaleManager(atc_->interscale_manager());
 
377
      const INT_ARRAY & nodeType((interscaleManager.dense_matrix_int("NodalGeometryType"))->quantity());
 
378
      if (atomicRegulator_->use_localized_lambda()) {
 
379
        for (int i = 0; i < nNodes_; ++i) {
 
380
          if (nodeType(i,0)==BOUNDARY) {
 
381
            for (int j = 0; j < nsd_; ++j) {
257
382
              hooverNodes_.insert(pair<int,int>(i,j));
 
383
            }
 
384
          }
 
385
        }
258
386
      }
259
387
      else {
260
 
        for (int i = 0; i < nNodes_; ++i)
261
 
          if (nodeType(i)==ATC_Transfer::BOUNDARY || nodeType(i)==ATC_Transfer::MD_ONLY)
262
 
            for (int j = 0; j < nsd_; ++j)
 
388
        for (int i = 0; i < nNodes_; ++i) {
 
389
          if (nodeType(i,0)==BOUNDARY || nodeType(i,0)==MD_ONLY) {
 
390
            for (int j = 0; j < nsd_; ++j) {
263
391
              hooverNodes_.insert(pair<int,int>(i,j));
 
392
            }
 
393
          }
 
394
        }
264
395
      }
265
396
    }
266
397
  }
270
401
  //            uses existing lambda to modify given
271
402
  //            atomic quantity
272
403
  //--------------------------------------------------------
273
 
  void GlcKinetostat::apply_to_atoms(double ** atomicQuantity,
 
404
  void GlcKinetostat::apply_to_atoms(PerAtomQuantity<double> * quantity,
274
405
                                     const DENS_MAT & lambdaAtom,
275
406
                                     double dt)
276
407
  {
277
 
    if (nLocal_>0) {
278
 
      for (int i = 0; i < nLocal_; i++) {
279
 
        for (int j = 0; j < nsd_; j++) {
280
 
          atomicQuantity[internalToAtom_(i)][j] -= lambdaAtom(i,j);
281
 
        }
282
 
      }
283
 
    }
284
 
  }
285
 
 
286
 
  //--------------------------------------------------------
287
 
  //  set_weights
288
 
  //            sets diagonal weighting matrix used in
289
 
  //            solve_for_lambda
290
 
  //--------------------------------------------------------
291
 
  void GlcKinetostat::set_weights(DIAG_MAT & weights)
292
 
  {
293
 
    if (nLocalLambda_>0) {
294
 
      DENS_VEC weightVector(nLocal_);
295
 
      weightVector = 1.;
296
 
      DENS_VEC maskedWeightVector(nLocalLambda_);
297
 
      maskedWeightVector = internalToOverlapMap_*weightVector;
298
 
      weights.reset(maskedWeightVector);
299
 
    }
300
 
  }
301
 
 
302
 
  //--------------------------------------------------------
303
 
  //  reset_nlocal:
304
 
  //    resets data dependent on local atom count
305
 
  //--------------------------------------------------------
306
 
  void GlcKinetostat::reset_nlocal()
307
 
  {
308
 
    RegulatorShapeFunction::reset_nlocal();
309
 
    if (nLocal_ > 0) {
310
 
      atcTransfer_->compute_atomic_mass(atomicMass_);
311
 
    }
312
 
  }
313
 
 
 
408
    *quantity -= lambdaAtom;
 
409
  }
 
410
 
314
411
  //--------------------------------------------------------
315
412
  //--------------------------------------------------------
316
413
  //  Class DisplacementGlc
323
420
  //--------------------------------------------------------
324
421
  DisplacementGlc::DisplacementGlc(Kinetostat * kinetostat) :
325
422
    GlcKinetostat(kinetostat),
326
 
    nodalDisplacements_(atcTransfer_->get_field(DISPLACEMENT)),
327
 
    x_(atcTransfer_->get_x())           
328
 
  {
 
423
    nodalAtomicMassWeightedDisplacement_(NULL),
 
424
    nodalDisplacements_(atc_->field(DISPLACEMENT))
 
425
  {
 
426
    // do nothing
 
427
  }
 
428
 
 
429
  //--------------------------------------------------------
 
430
  //  constructor_transfers
 
431
  //    instantiates or obtains all dependency managed data
 
432
  //--------------------------------------------------------
 
433
  void DisplacementGlc::construct_transfers()
 
434
  {
 
435
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
436
 
 
437
    // set up node mappings
 
438
    create_node_maps();
 
439
 
 
440
    // set up shape function matrix
 
441
    if (this->use_local_shape_functions()) {
 
442
      lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
 
443
      interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
 
444
                                                  regulatorPrefix_+"LambdaAtomMap");
 
445
      shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
 
446
                                                           lambdaAtomMap_,
 
447
                                                           nodeToOverlapMap_);
 
448
    }
 
449
    else {
 
450
      shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
 
451
    }
 
452
    interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
 
453
                                                 regulatorPrefix_+"LambdaCouplingMatrixMomentum");
 
454
 
 
455
    // set linear solver strategy
 
456
    if (atomicRegulator_->use_lumped_lambda_solve()) {
 
457
      linearSolverType_ = AtomicRegulator::RSL_SOLVE;
 
458
    }
 
459
    else {
 
460
      linearSolverType_ = AtomicRegulator::CG_SOLVE;
 
461
    }
 
462
    
 
463
 
 
464
    // base class transfers
 
465
    GlcKinetostat::construct_transfers();
 
466
 
 
467
    // atomic force induced by kinetostat
 
468
    atomKinetostatForce_ = new AtomicKinetostatForceDisplacement(atc_);
 
469
    interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
 
470
                                            regulatorPrefix_+"AtomKinetostatForce");
 
471
        
 
472
    // restricted force due to kinetostat
 
473
    nodalAtomicLambdaForce_ = new AtfShapeFunctionRestriction(atc_,
 
474
                                                              atomKinetostatForce_,
 
475
                                                              interscaleManager.per_atom_sparse_matrix("Interpolant"));
 
476
    interscaleManager.add_dense_matrix(nodalAtomicLambdaForce_,
 
477
                                       regulatorPrefix_+"NodalAtomicLambdaForce");
 
478
    
 
479
    // nodal displacement restricted from atoms
 
480
    nodalAtomicMassWeightedDisplacement_ = interscaleManager.dense_matrix("NodalAtomicMassWeightedDisplacement");
 
481
  }
 
482
 
 
483
  //--------------------------------------------------------
 
484
  //  initialize
 
485
  //    initializes all method data
 
486
  //--------------------------------------------------------
 
487
  void DisplacementGlc::initialize()
 
488
  {
 
489
    GlcKinetostat::initialize();
 
490
 
329
491
    // sets up time filter for cases where variables temporally filtered
330
 
    TimeFilterManager * timeFilterManager = atcTransfer_->get_time_filter_manager();
 
492
    TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
331
493
    if (!timeFilterManager->end_equilibrate()) {
332
 
      nodalAtomicLambdaForce_ = 0.;
333
 
      lambdaForceFiltered_ = 0.;
334
 
      timeFilter_->initialize(lambdaForceFiltered_);
 
494
      *lambdaForceFiltered_ = 0.;
 
495
      timeFilter_->initialize(lambdaForceFiltered_->quantity());
335
496
    }
336
497
  }
337
 
 
 
498
 
338
499
  //--------------------------------------------------------
339
500
  //  apply:
340
501
  //    apply the kinetostat to the atoms
359
520
    set_kinetostat_rhs(rhs,dt);
360
521
 
361
522
    // solve linear system for lambda
362
 
    solve_for_lambda(rhs);
363
 
 
364
 
    // compute force applied by lambda
365
 
    DENS_MAT lambdaAtom;
366
 
    compute_lambda_force(lambdaForce_,lambdaAtom,dt);
 
523
    solve_for_lambda(rhs,lambda_->set_quantity());
367
524
 
368
525
    // compute nodal atomic power
369
526
    compute_nodal_lambda_force(dt);
370
527
 
371
528
    // apply kinetostat to atoms
372
 
    apply_to_atoms(x_,lambdaAtom);
 
529
    apply_to_atoms(atomPositions_,atomLambdas_->quantity());
373
530
  }
374
531
 
375
532
  //--------------------------------------------------------
380
537
  void DisplacementGlc::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
381
538
  {
382
539
    // form rhs : sum_a (hatN_Ia * x_ai) - (Upsilon)_Ii
383
 
    DENS_MAT atomicDisplacement;
384
 
    atcTransfer_->compute_atomic_centerOfMass_displacement(atomicDisplacement,x_);
385
 
    rhs.reset(nNodes_,nsd_);
386
 
    atcTransfer_->restrict_volumetric_quantity(atomicDisplacement,rhs);
387
 
    rhs -= (atcTransfer_->get_mass_mat_md(VELOCITY))*nodalDisplacements_;
388
 
  }
389
 
 
390
 
  //--------------------------------------------------------
391
 
  //  compute_lambda_force
392
 
  //            compute the equivalent force on the atoms
393
 
  //            induced by lambda
394
 
  //--------------------------------------------------------
395
 
  void DisplacementGlc::compute_lambda_force(DENS_MAT & lambdaForce,
396
 
                                             DENS_MAT & lambdaAtom,
397
 
                                             double dt)
398
 
  {
399
 
    if (nLocal_>0) {
400
 
      // prolongation to (unique) nodes
401
 
      lambdaAtom.reset(nLocal_,nsd_);
402
 
      atcTransfer_->prolong(lambda_,lambdaAtom);
403
 
 
404
 
      for (int i = 0; i < nLocal_; i++)
405
 
        for (int j = 0; j < nsd_; j++)
406
 
          lambdaForce(i,j) =  (-1./(dt*dt))*lambdaAtom(i,j)*atomicMass_(i);
407
 
    }
 
540
    rhs = nodalAtomicMassWeightedDisplacement_->quantity();
 
541
    rhs -= ((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalDisplacements_.quantity());
408
542
  }
409
543
 
410
544
  //--------------------------------------------------------
414
548
  //--------------------------------------------------------
415
549
  void DisplacementGlc::compute_nodal_lambda_force(double dt)
416
550
  {
417
 
    atcTransfer_->restrict_volumetric_quantity(lambdaForce_,nodalAtomicLambdaForce_);
418
 
    timeFilter_->apply_post_step1(lambdaForceFiltered_,nodalAtomicLambdaForce_,dt);
 
551
    const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
 
552
    timeFilter_->apply_post_step1(lambdaForceFiltered_->set_quantity(),
 
553
                                  myNodalAtomicLambdaForce,dt);
419
554
 
420
555
    // update FE displacements for localized thermostats
421
 
    apply_localization_correction(nodalAtomicLambdaForce_,
422
 
                                  nodalDisplacements_,
 
556
    apply_localization_correction(myNodalAtomicLambdaForce,
 
557
                                  nodalDisplacements_.set_quantity(),
423
558
                                  dt*dt);
424
559
  }
425
560
 
432
567
  {
433
568
    // apply time filtered lambda force
434
569
    DENS_MAT lambdaZero(nNodes_,nsd_);
435
 
    timeFilter_->apply_pre_step1(lambdaForceFiltered_,(-1./dt/dt)*lambdaZero,dt);
 
570
    timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),(-1./dt/dt)*lambdaZero,dt);
436
571
  }
437
572
 
438
573
  //--------------------------------------------------------
440
575
  //            sets diagonal weighting matrix used in
441
576
  //            solve_for_lambda
442
577
  //--------------------------------------------------------
443
 
  void DisplacementGlc::set_weights(DIAG_MAT & weights)
 
578
  void DisplacementGlc::set_weights()
444
579
  {
445
 
    if (nLocalLambda_>0) {
446
 
      DENS_VEC maskedWeightVector(nLocalLambda_);
447
 
      maskedWeightVector = internalToOverlapMap_*atomicMass_;
448
 
      weights.reset(maskedWeightVector);
 
580
    if (lambdaAtomMap_) {
 
581
      MappedAtomQuantity * myWeights = new MappedAtomQuantity(atc_,atomMasses_,lambdaAtomMap_);
 
582
      weights_ = myWeights;
 
583
      (atc_->interscale_manager()).add_per_atom_quantity(myWeights,
 
584
                                                         "AtomMassesMapped");
 
585
    }
 
586
    else {
 
587
      weights_ = atomMasses_;
449
588
    }
450
589
  }
451
590
 
459
598
                                                      DENS_MAT & nodalField,
460
599
                                                      double weight)
461
600
  {
462
 
    // NOTE can add check to see if set is empty for faster performance
 
601
    
463
602
    DENS_MAT nodalLambdaRoc(nNodes_,nsd_);
464
 
    atcTransfer_->apply_inverse_mass_matrix(source,
465
 
                                            nodalLambdaRoc,
466
 
                                            VELOCITY);
 
603
    atc_->apply_inverse_mass_matrix(source,
 
604
                                    nodalLambdaRoc,
 
605
                                    VELOCITY);
467
606
    set<pair<int,int> >::const_iterator iter;
468
607
    for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
469
608
      nodalLambdaRoc(iter->first,iter->second) = 0.;
476
615
  //  output:
477
616
  //    adds all relevant output to outputData
478
617
  //--------------------------------------------------------
479
 
  void DisplacementGlc::output(double dt, OUTPUT_LIST & outputData)
 
618
  void DisplacementGlc::output(OUTPUT_LIST & outputData)
480
619
  {
481
 
    outputData["lambda"] = &(lambda_);
482
 
    outputData["nodalLambdaForce"] = &(nodalAtomicLambdaForce_);
 
620
    _nodalAtomicLambdaForceOut_ = nodalAtomicLambdaForce_->quantity();
 
621
    DENS_MAT & lambda(lambda_->set_quantity());
 
622
    if ((atc_->lammps_interface())->rank_zero()) {
 
623
      outputData[regulatorPrefix_+"Lambda"] = &lambda;
 
624
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_nodalAtomicLambdaForceOut_);
 
625
    }
483
626
  }
484
627
  
485
628
  //--------------------------------------------------------
494
637
  //--------------------------------------------------------
495
638
  DisplacementGlcFiltered::DisplacementGlcFiltered(Kinetostat * kinetostat) :
496
639
    DisplacementGlc(kinetostat),
497
 
    nodalAtomicDisplacements_(atcTransfer_->get_atomic_field(DISPLACEMENT))
 
640
    nodalAtomicDisplacements_(atc_->nodal_atomic_field(DISPLACEMENT))
498
641
  {
499
642
    // do nothing
500
643
  }
509
652
    // apply time filtered lambda to atomic fields
510
653
    DisplacementGlc::apply_pre_filtering(dt);
511
654
    DENS_MAT nodalAcceleration(nNodes_,nsd_);
512
 
    atcTransfer_->apply_inverse_md_mass_matrix(lambdaForceFiltered_,
513
 
                                               nodalAcceleration,
514
 
                                               VELOCITY);
 
655
    atc_->apply_inverse_md_mass_matrix(lambdaForceFiltered_->set_quantity(),
 
656
                                       nodalAcceleration,
 
657
                                       VELOCITY);
515
658
    nodalAtomicDisplacements_ += dt*dt*nodalAcceleration;
516
659
  }
517
660
 
523
666
  void DisplacementGlcFiltered::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
524
667
  {
525
668
    // form rhs : sum_a (hatN_Ia * x_ai) - (Upsilon)_Ii
526
 
    double coef = 1./(timeFilter_->get_unfiltered_coefficient_pre_s1(dt));
527
 
    rhs = coef*(atcTransfer_->get_mass_mat_md(VELOCITY))*(nodalAtomicDisplacements_ - nodalDisplacements_);
 
669
    double coef = 1./(timeFilter_->unfiltered_coefficient_pre_s1(dt));
 
670
    rhs = coef*((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalAtomicDisplacements_.quantity() - nodalDisplacements_.quantity());
528
671
  }
529
672
 
530
673
  //--------------------------------------------------------
534
677
  //--------------------------------------------------------
535
678
  void DisplacementGlcFiltered::compute_nodal_lambda_force(double dt)
536
679
  {
537
 
    atcTransfer_->restrict_volumetric_quantity(lambdaForce_,nodalAtomicLambdaForce_);
538
 
    timeFilter_->apply_post_step1(lambdaForceFiltered_,nodalAtomicLambdaForce_,dt);
 
680
    const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
 
681
    DENS_MAT & myLambdaForceFiltered(lambdaForceFiltered_->set_quantity());
 
682
    timeFilter_->apply_post_step1(myLambdaForceFiltered,
 
683
                                  myNodalAtomicLambdaForce,dt);
539
684
 
540
685
    // update filtered atomic displacements
541
 
    DENS_MAT nodalLambdaRoc(nodalAtomicLambdaForce_.nRows(),nodalAtomicLambdaForce_.nCols());
542
 
    atcTransfer_->apply_inverse_md_mass_matrix(nodalAtomicLambdaForce_,
543
 
                                               nodalLambdaRoc,
544
 
                                               VELOCITY);
545
 
    timeFilter_->apply_post_step1(nodalAtomicDisplacements_,dt*dt*nodalLambdaRoc,dt);
 
686
    DENS_MAT nodalLambdaRoc(myNodalAtomicLambdaForce.nRows(),myNodalAtomicLambdaForce.nCols());
 
687
    atc_->apply_inverse_md_mass_matrix(myNodalAtomicLambdaForce,
 
688
                                       nodalLambdaRoc,
 
689
                                       VELOCITY);
 
690
    timeFilter_->apply_post_step1(nodalAtomicDisplacements_.set_quantity(),dt*dt*nodalLambdaRoc,dt);
546
691
 
547
692
    // update FE displacements for localized thermostats
548
 
    apply_localization_correction(lambdaForceFiltered_,
549
 
                                  nodalDisplacements_,
 
693
    apply_localization_correction(myLambdaForceFiltered,
 
694
                                  nodalDisplacements_.set_quantity(),
550
695
                                  dt*dt);
551
696
  }
552
697
 
554
699
  //  output:
555
700
  //    adds all relevant output to outputData
556
701
  //--------------------------------------------------------
557
 
  void DisplacementGlcFiltered::output(double dt, OUTPUT_LIST & outputData)
 
702
  void DisplacementGlcFiltered::output(OUTPUT_LIST & outputData)
558
703
  {
559
 
    outputData["lambda"] = &(lambda_);
560
 
    outputData["nodalLambdaForce"] = &(lambdaForceFiltered_);
 
704
    DENS_MAT & lambda(lambda_->set_quantity());
 
705
    DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
 
706
    if ((atc_->lammps_interface())->rank_zero()) {
 
707
      outputData[regulatorPrefix_+"Lambda"] = &lambda;
 
708
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &lambdaForceFiltered;
 
709
    }
561
710
  }
562
711
  
563
712
  //--------------------------------------------------------
572
721
  //--------------------------------------------------------
573
722
  VelocityGlc::VelocityGlc(Kinetostat * kinetostat) :
574
723
    GlcKinetostat(kinetostat),
575
 
    nodalVelocities_(atcTransfer_->get_field(VELOCITY)),
576
 
    v_(atcTransfer_->get_v())
577
 
  {
 
724
    nodalAtomicMomentum_(NULL),
 
725
    nodalVelocities_(atc_->field(VELOCITY))
 
726
  {
 
727
    // do nothing
 
728
  }
 
729
 
 
730
  //--------------------------------------------------------
 
731
  //  constructor_transfers
 
732
  //    instantiates or obtains all dependency managed data
 
733
  //--------------------------------------------------------
 
734
  void VelocityGlc::construct_transfers()
 
735
  {
 
736
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
737
 
 
738
    // set up node mappings
 
739
    create_node_maps();
 
740
 
 
741
    // set up shape function matrix
 
742
    if (this->use_local_shape_functions()) {
 
743
      lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
 
744
      interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
 
745
                                                  regulatorPrefix_+"LambdaAtomMap");
 
746
      shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
 
747
                                                           lambdaAtomMap_,
 
748
                                                           nodeToOverlapMap_);
 
749
    }
 
750
    else {
 
751
      shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
 
752
    }
 
753
    interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
 
754
                                                 regulatorPrefix_+"LambdaCouplingMatrixMomentum");
 
755
 
 
756
    // set linear solver strategy
 
757
    if (atomicRegulator_->use_lumped_lambda_solve()) {
 
758
      linearSolverType_ = AtomicRegulator::RSL_SOLVE;
 
759
    }
 
760
    else {
 
761
      linearSolverType_ = AtomicRegulator::CG_SOLVE;
 
762
    }
 
763
    
 
764
    // base class transfers
 
765
    GlcKinetostat::construct_transfers();
 
766
 
 
767
    // atomic force induced by kinetostat
 
768
    atomKinetostatForce_ = new AtomicKinetostatForceVelocity(atc_);
 
769
    interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
 
770
                                            regulatorPrefix_+"AtomKinetostatForce");
 
771
        
 
772
    // restricted force due to kinetostat
 
773
    nodalAtomicLambdaForce_ = new AtfShapeFunctionRestriction(atc_,
 
774
                                                              atomKinetostatForce_,
 
775
                                                              interscaleManager.per_atom_sparse_matrix("Interpolant"));
 
776
    interscaleManager.add_dense_matrix(nodalAtomicLambdaForce_,
 
777
                                            regulatorPrefix_+"NodalAtomicLambdaForce");
 
778
    
 
779
    // nodal momentum restricted from atoms
 
780
    nodalAtomicMomentum_ = interscaleManager.dense_matrix("NodalAtomicMomentum");
 
781
 
 
782
    
 
783
  }
 
784
 
 
785
  //--------------------------------------------------------
 
786
  //  initialize
 
787
  //    initializes all method data
 
788
  //--------------------------------------------------------
 
789
  void VelocityGlc::initialize()
 
790
  {
 
791
    GlcKinetostat::initialize();
 
792
 
578
793
    // sets up time filter for cases where variables temporally filtered
579
 
    TimeFilterManager * timeFilterManager = atcTransfer_->get_time_filter_manager();
 
794
    TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
580
795
    if (!timeFilterManager->end_equilibrate()) {
581
 
      nodalAtomicLambdaForce_ = 0.;
582
 
      lambdaForceFiltered_ = 0.;
583
 
      timeFilter_->initialize(lambdaForceFiltered_);
 
796
      lambdaForceFiltered_->set_quantity() = 0.;
 
797
      timeFilter_->initialize(lambdaForceFiltered_->quantity());
584
798
    }
585
799
  }
586
 
 
 
800
 
587
801
  //--------------------------------------------------------
588
802
  //  apply_mid_corrector:
589
803
  //    apply the kinetostat during the middle of the
602
816
  void VelocityGlc::apply_post_corrector(double dt)
603
817
  {
604
818
    double dtLambda = 0.5*dt;
605
 
    compute_kinetostat(dt);
 
819
    compute_kinetostat(dtLambda);
606
820
  }
607
821
 
608
822
  //--------------------------------------------------------
614
828
  {
615
829
    // apply time filtered lambda to atomic fields
616
830
    DENS_MAT lambdaZero(nNodes_,nsd_);
617
 
    timeFilter_->apply_pre_step1(lambdaForceFiltered_,(-1./dt)*lambdaZero,dt);
 
831
    timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),(-1./dt)*lambdaZero,dt);
618
832
  }
619
833
 
620
834
  //--------------------------------------------------------
632
846
    set_kinetostat_rhs(rhs,dt);
633
847
    
634
848
    // solve linear system for lambda
635
 
    solve_for_lambda(rhs);
636
 
    
637
 
    // compute force applied by lambda
638
 
    DENS_MAT lambdaAtom;
639
 
    compute_lambda_force(lambdaForce_,lambdaAtom,dt);
 
849
    solve_for_lambda(rhs,lambda_->set_quantity());
640
850
 
641
851
    // compute nodal atomic power
642
852
    compute_nodal_lambda_force(dt);
643
853
 
644
854
    // apply kinetostat to atoms
645
 
    apply_to_atoms(v_,lambdaAtom);
 
855
    apply_to_atoms(atomVelocities_,atomLambdas_->quantity());
646
856
  }
647
857
 
648
858
  //--------------------------------------------------------
653
863
  void VelocityGlc::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
654
864
  {
655
865
    // form rhs : sum_a (hatN_Ia * x_ai) - (\dot{Upsilon})_Ii
656
 
    DENS_MAT atomicMomentum;
657
 
    atcTransfer_->compute_atomic_momentum(atomicMomentum,v_);
658
 
    rhs.reset(nNodes_,nsd_);
659
 
    atcTransfer_->restrict_volumetric_quantity(atomicMomentum,rhs);
660
 
    rhs -= (atcTransfer_->get_mass_mat_md(VELOCITY))*nodalVelocities_;
661
 
  }
662
 
 
663
 
  //--------------------------------------------------------
664
 
  //  compute_lambda_force
665
 
  //            compute the equivalent force on the atoms
666
 
  //            induced by lambda
667
 
  //--------------------------------------------------------
668
 
  void VelocityGlc::compute_lambda_force(DENS_MAT & lambdaForce,
669
 
                                          DENS_MAT & lambdaAtom,
670
 
                                          double dt)
671
 
  {
672
 
    if (nLocal_>0) {
673
 
      // prolongation to (unique) nodes
674
 
      lambdaAtom.reset(nLocal_,nsd_);
675
 
      atcTransfer_->prolong(lambda_,lambdaAtom);
676
 
 
677
 
      for (int i = 0; i < nLocal_; i++)
678
 
        for (int j = 0; j < nsd_; j++)
679
 
          lambdaForce(i,j) =  (-1./(dt))*lambdaAtom(i,j)*atomicMass_(i);
680
 
    }
 
866
    rhs = nodalAtomicMomentum_->quantity();
 
867
    rhs -= ((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalVelocities_.quantity());
681
868
  }
682
869
 
683
870
  //--------------------------------------------------------
687
874
  //--------------------------------------------------------
688
875
  void VelocityGlc::compute_nodal_lambda_force(double dt)
689
876
  {
690
 
    atcTransfer_->restrict_volumetric_quantity(lambdaForce_,nodalAtomicLambdaForce_);
691
 
    timeFilter_->apply_pre_step1(lambdaForceFiltered_,nodalAtomicLambdaForce_,dt);
 
877
    const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
 
878
 
 
879
    timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),
 
880
                                 myNodalAtomicLambdaForce,dt);
692
881
 
693
882
    // update FE displacements for localized thermostats
694
 
    apply_localization_correction(nodalAtomicLambdaForce_,
695
 
                                  nodalVelocities_,
 
883
    apply_localization_correction(myNodalAtomicLambdaForce,
 
884
                                  nodalVelocities_.set_quantity(),
696
885
                                  dt);
697
886
  }
698
887
 
701
890
  //            sets diagonal weighting matrix used in
702
891
  //            solve_for_lambda
703
892
  //--------------------------------------------------------
704
 
  void VelocityGlc::set_weights(DIAG_MAT & weights)
 
893
  void VelocityGlc::set_weights()
705
894
  {
706
 
    if (nLocalLambda_>0) {
707
 
      DENS_VEC maskedWeightVector = internalToOverlapMap_*atomicMass_;
708
 
      weights.reset(maskedWeightVector);
 
895
    if (lambdaAtomMap_) {
 
896
      MappedAtomQuantity * myWeights = new MappedAtomQuantity(atc_,atomMasses_,lambdaAtomMap_);
 
897
      weights_ = myWeights;
 
898
      (atc_->interscale_manager()).add_per_atom_quantity(myWeights,
 
899
                                                         "AtomMassesMapped");
 
900
    }
 
901
    else {
 
902
      weights_ = atomMasses_;
709
903
    }
710
904
  }
711
905
 
719
913
                                                  DENS_MAT & nodalField,
720
914
                                                  double weight)
721
915
  {
722
 
    // NOTE can add check to see if set is empty for faster performance
 
916
    
723
917
    DENS_MAT nodalLambdaRoc(nNodes_,nsd_);
724
 
    atcTransfer_->apply_inverse_mass_matrix(source,
725
 
                                            nodalLambdaRoc,
726
 
                                            VELOCITY);
 
918
    atc_->apply_inverse_mass_matrix(source,
 
919
                                    nodalLambdaRoc,
 
920
                                    VELOCITY);
727
921
    set<pair<int,int> >::const_iterator iter;
728
922
    for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
729
923
      nodalLambdaRoc(iter->first,iter->second) = 0.;
730
924
    }
731
 
 
 
925
   
732
926
    nodalField += weight*nodalLambdaRoc;
733
927
  }
734
928
 
736
930
  //  output:
737
931
  //    adds all relevant output to outputData
738
932
  //--------------------------------------------------------
739
 
  void VelocityGlc::output(double dt, OUTPUT_LIST & outputData)
 
933
  void VelocityGlc::output(OUTPUT_LIST & outputData)
740
934
  {
741
 
    outputData["lambda"] = &(lambda_);
742
 
    outputData["nodalLambdaForce"] = &(nodalAtomicLambdaForce_);
 
935
    _nodalAtomicLambdaForceOut_ = nodalAtomicLambdaForce_->quantity();
 
936
    if ((atc_->lammps_interface())->rank_zero()) {
 
937
      outputData[regulatorPrefix_+"Lambda"] = &(lambda_->set_quantity());
 
938
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_nodalAtomicLambdaForceOut_);
 
939
    }
743
940
  }
744
941
 
745
942
  //--------------------------------------------------------
754
951
  //--------------------------------------------------------
755
952
  VelocityGlcFiltered::VelocityGlcFiltered(Kinetostat *kinetostat) 
756
953
    : VelocityGlc(kinetostat),
757
 
      nodalAtomicVelocities_(atcTransfer_->get_atomic_field(VELOCITY))
 
954
      nodalAtomicVelocities_(atc_->nodal_atomic_field(VELOCITY))
758
955
  {
759
956
    // do nothing
760
957
  }
769
966
    // apply time filtered lambda to atomic fields
770
967
    VelocityGlc::apply_pre_filtering(dt);
771
968
    DENS_MAT nodalAcceleration(nNodes_,nsd_);
772
 
    atcTransfer_->apply_inverse_md_mass_matrix(lambdaForceFiltered_,
773
 
                                               nodalAcceleration,
774
 
                                               VELOCITY);
 
969
    atc_->apply_inverse_md_mass_matrix(lambdaForceFiltered_->quantity(),
 
970
                                       nodalAcceleration,
 
971
                                       VELOCITY);
775
972
    nodalAtomicVelocities_ += dt*nodalAcceleration;
776
973
  }
777
974
 
783
980
  void VelocityGlcFiltered::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
784
981
  {
785
982
    // form rhs : sum_a (hatN_Ia * x_ai) - (Upsilon)_Ii
786
 
    double coef = 1./(timeFilter_->get_unfiltered_coefficient_pre_s1(dt));
787
 
    rhs = coef*(atcTransfer_->get_mass_mat_md(VELOCITY))*(nodalAtomicVelocities_ - nodalVelocities_);
 
983
    double coef = 1./(timeFilter_->unfiltered_coefficient_pre_s1(dt));
 
984
    rhs = coef*((atc_->mass_mat_md(VELOCITY)).quantity())*(nodalAtomicVelocities_.quantity() - nodalVelocities_.quantity());
788
985
  }
789
986
 
790
987
  //--------------------------------------------------------
794
991
  //--------------------------------------------------------
795
992
  void VelocityGlcFiltered::compute_nodal_lambda_force(double dt)
796
993
  {
797
 
    atcTransfer_->restrict_volumetric_quantity(lambdaForce_,nodalAtomicLambdaForce_);
798
 
    timeFilter_->apply_post_step1(lambdaForceFiltered_,nodalAtomicLambdaForce_,dt);
 
994
    const DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->quantity());
 
995
    DENS_MAT & myLambdaForceFiltered(lambdaForceFiltered_->set_quantity());
 
996
    timeFilter_->apply_post_step1(myLambdaForceFiltered,myNodalAtomicLambdaForce,dt);
799
997
 
800
998
    // update filtered atomic displacements
801
 
    DENS_MAT nodalLambdaRoc(nodalAtomicLambdaForce_.nRows(),nodalAtomicLambdaForce_.nCols());
802
 
    atcTransfer_->apply_inverse_md_mass_matrix(nodalAtomicLambdaForce_,
803
 
                                               nodalLambdaRoc,
804
 
                                               VELOCITY);
805
 
    timeFilter_->apply_post_step1(nodalAtomicVelocities_,dt*nodalLambdaRoc,dt);
 
999
    DENS_MAT nodalLambdaRoc(myNodalAtomicLambdaForce.nRows(),myNodalAtomicLambdaForce.nCols());
 
1000
    atc_->apply_inverse_md_mass_matrix(myNodalAtomicLambdaForce,
 
1001
                                       nodalLambdaRoc,
 
1002
                                       VELOCITY);
 
1003
    timeFilter_->apply_post_step1(nodalAtomicVelocities_.set_quantity(),dt*nodalLambdaRoc,dt);
806
1004
 
807
1005
    // update FE displacements for localized thermostats
808
 
    apply_localization_correction(lambdaForceFiltered_,
809
 
                                 nodalVelocities_,
810
 
                                 dt);
 
1006
    apply_localization_correction(myLambdaForceFiltered,
 
1007
                                  nodalVelocities_.set_quantity(),
 
1008
                                  dt);
811
1009
  }
812
1010
 
813
1011
  //--------------------------------------------------------
814
1012
  //  output:
815
1013
  //    adds all relevant output to outputData
816
1014
  //--------------------------------------------------------
817
 
  void VelocityGlcFiltered::output(double dt, OUTPUT_LIST & outputData)
 
1015
  void VelocityGlcFiltered::output(OUTPUT_LIST & outputData)
818
1016
  {
819
 
    outputData["lambda"] = &(lambda_);
820
 
    outputData["nodalLambdaForce"] = &(lambdaForceFiltered_);
 
1017
    if ((atc_->lammps_interface())->rank_zero()) {
 
1018
      outputData[regulatorPrefix_+"Lambda"] = &(lambda_->set_quantity());
 
1019
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &(lambdaForceFiltered_->set_quantity());
 
1020
    }
821
1021
  }
822
1022
 
823
1023
  //--------------------------------------------------------
832
1032
  //--------------------------------------------------------
833
1033
  StressFlux::StressFlux(Kinetostat * kinetostat) :
834
1034
    GlcKinetostat(kinetostat),
835
 
    nodalForce_(atcTransfer_->get_field_rhs(VELOCITY)),
836
 
    nodalAtomicForce_(atcTransfer_->get_fe_atomic_field_roc(VELOCITY)),
837
 
    momentumSource_(atcTransfer_->get_atomic_source(VELOCITY)),
838
 
    v_(atcTransfer_->get_v()),
839
 
    f_(atcTransfer_->get_f())
 
1035
    nodalForce_(atc_->field_rhs(VELOCITY)),
 
1036
    nodalAtomicForce_(NULL),
 
1037
    nodalGhostForce_(NULL),
 
1038
    momentumSource_(atc_->atomic_source(VELOCITY))
840
1039
  {
841
1040
    // flag for performing boundary flux calculation
842
 
    if (kinetostat_->get_coupling_mode()==Kinetostat::FLUX)
843
 
      fieldMask_(VELOCITY,FLUX) = true;
844
 
 
845
 
    // sets up space for ghost force related variables
846
 
    nodalGhostForce_.reset(nNodes_,nsd_);
847
 
    nodalGhostForceFiltered_.reset(nNodes_,nsd_);
848
 
 
849
 
    // NOTE ifdefs are for using reference lammps force as base pressure
850
 
#if false
851
 
    int nLocalTotal = atcTransfer->get_nlocal_total();
852
 
    int nsd = atcTransfer->get_nsd();
853
 
    f0_ = new double*[nLocalTotal];
854
 
    for (int i = 0; i < nLocalTotal; ++i) {
855
 
      f0_[i] = new double[nsd];
856
 
      for (int j = 0; j < nsd; ++j)
857
 
        f0_[i][j] = f_[i][j];
858
 
    }
859
 
#endif
 
1041
    fieldMask_(VELOCITY,FLUX) = true;
860
1042
  }
861
1043
 
862
1044
  StressFlux::~StressFlux()
863
1045
  {
864
 
#if false
865
 
    if (f0_) {
866
 
      ATC_Transfer * atcTransfer = kinetostat_->get_atc_transfer();
867
 
      int nLocalTotal = atcTransfer->get_nlocal_total();
868
 
      for (int i = 0; i < nLocalTotal; ++i)
869
 
        delete [] f0_[i];
870
 
      delete [] f0_;
871
 
    }
872
 
#endif
 
1046
    // do nothing
 
1047
  }
 
1048
 
 
1049
  //--------------------------------------------------------
 
1050
  //  constructor_transfers
 
1051
  //    instantiates or obtains all dependency managed data
 
1052
  //--------------------------------------------------------
 
1053
  void StressFlux::construct_transfers()
 
1054
  {
 
1055
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
1056
 
 
1057
    // set up node mappings
 
1058
    create_node_maps();
 
1059
 
 
1060
    // set up shape function matrix
 
1061
    if (this->use_local_shape_functions()) {
 
1062
      lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
 
1063
      interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
 
1064
                                                  regulatorPrefix_+"LambdaAtomMap");
 
1065
      shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
 
1066
                                                           lambdaAtomMap_,
 
1067
                                                           nodeToOverlapMap_);
 
1068
    }
 
1069
    else {
 
1070
      shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
 
1071
    }
 
1072
    interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
 
1073
                                                 regulatorPrefix_+"LambdaCouplingMatrixMomentum");
 
1074
 
 
1075
    // set linear solver strategy
 
1076
    if (atomicRegulator_->use_lumped_lambda_solve()) {
 
1077
      linearSolverType_ = AtomicRegulator::RSL_SOLVE;
 
1078
    }
 
1079
    else {
 
1080
      linearSolverType_ = AtomicRegulator::CG_SOLVE;
 
1081
    }
 
1082
 
 
1083
    // base class transfers
 
1084
    GlcKinetostat::construct_transfers();
 
1085
 
 
1086
    // force at nodes due to atoms
 
1087
    nodalAtomicForce_ = interscaleManager.dense_matrix("NodalAtomicForce");
 
1088
 
 
1089
    // atomic force induced by kinetostat
 
1090
    atomKinetostatForce_ = new AtomicKinetostatForceStress(atc_,atomLambdas_);
 
1091
    interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
 
1092
                                            regulatorPrefix_+"AtomKinetostatForce");
 
1093
        
 
1094
    // restricted force due to kinetostat
 
1095
    nodalAtomicLambdaForce_ = new AtfShapeFunctionRestriction(atc_,
 
1096
                                                              atomKinetostatForce_,
 
1097
                                                              interscaleManager.per_atom_sparse_matrix("Interpolant"));
 
1098
    interscaleManager.add_dense_matrix(nodalAtomicLambdaForce_,
 
1099
                                       regulatorPrefix_+"NodalAtomicLambdaForce");
 
1100
 
 
1101
    // sets up space for ghost force related variables
 
1102
    if (atc_->groupbit_ghost()) {
 
1103
      GhostCouplingMatrix * shapeFunctionGhost = new GhostCouplingMatrix(atc_,interscaleManager.per_atom_sparse_matrix("InterpolantGhost"),
 
1104
                                                                         regulatedNodes_,nodeToOverlapMap_);
 
1105
      interscaleManager.add_sparse_matrix(shapeFunctionGhost,
 
1106
                                          regulatorPrefix_+"GhostCouplingMatrix");
 
1107
      FundamentalAtomQuantity * atomGhostForce = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE,
 
1108
                                                                                                 GHOST);
 
1109
      nodalGhostForce_ = new AtfShapeFunctionRestriction(atc_,atomGhostForce,
 
1110
                                                         shapeFunctionGhost);
 
1111
      interscaleManager.add_dense_matrix(nodalGhostForce_,
 
1112
                                         regulatorPrefix_+"NodalGhostForce");
 
1113
      nodalGhostForceFiltered_.reset(nNodes_,nsd_);
 
1114
    }
 
1115
  }
 
1116
 
 
1117
  //--------------------------------------------------------
 
1118
  //  compute_boundary_flux:
 
1119
  //    computes the boundary flux to be consistent with
 
1120
  //    the controller
 
1121
  //--------------------------------------------------------
 
1122
  void StressFlux::compute_boundary_flux(FIELDS & fields)
 
1123
  {
 
1124
    GlcKinetostat::compute_boundary_flux(fields);
873
1125
  }
874
1126
 
875
1127
  //--------------------------------------------------------
877
1129
  //    apply the kinetostat to the atoms in the
878
1130
  //    mid-predictor integration phase
879
1131
  //--------------------------------------------------------
880
 
  void StressFlux::apply_mid_predictor(double dt)
 
1132
  void StressFlux::apply_pre_predictor(double dt)
881
1133
  {
882
1134
    double dtLambda = 0.5*dt;
883
1135
    // apply lambda force to atoms
884
 
    apply_to_atoms(v_,lambdaForce_,dtLambda);
 
1136
    apply_to_atoms(atomVelocities_,atomKinetostatForce_->quantity(),dtLambda);
885
1137
  }
886
1138
 
887
1139
  //--------------------------------------------------------
888
 
  //  apply_pre_corrector:
889
 
  //    apply the kinetostat to the atoms in the
890
 
  //    pre-corrector integration phase
891
 
  //--------------------------------------------------------
892
 
  void StressFlux::apply_pre_corrector(double dt)
893
 
  {
894
 
    // compute the kinetostat force
895
 
    compute_kinetostat(dt);
896
 
  }
897
 
 
898
 
  //--------------------------------------------------------
899
1140
  //  apply_post_corrector:
900
1141
  //    apply the kinetostat to the atoms in the
901
1142
  //    post-corrector integration phase
904
1145
  {
905
1146
    double dtLambda = 0.5*dt;
906
1147
    // apply lambda force to atoms
907
 
    apply_to_atoms(v_,lambdaForce_,dtLambda);
 
1148
    apply_to_atoms(atomVelocities_,atomKinetostatForce_->quantity(),dtLambda);
908
1149
  }
909
1150
 
910
1151
  //--------------------------------------------------------
914
1155
  //--------------------------------------------------------
915
1156
  void StressFlux::compute_kinetostat(double dt)
916
1157
  {
917
 
    // set up ghost force
918
 
    compute_ghost_force(nodalGhostForce_);
919
 
 
920
1158
    // initial filtering update
921
1159
    apply_pre_filtering(dt);
922
1160
 
923
1161
    // set up rhs
924
1162
    DENS_MAT rhs(nNodes_,nsd_);
925
1163
    set_kinetostat_rhs(rhs,dt);
926
 
    
 
1164
 
927
1165
    // solve linear system for lambda
928
 
    solve_for_lambda(rhs);
929
 
    
930
 
    // compute force applied by lambda
931
 
    compute_lambda_force(lambdaForce_);
 
1166
    solve_for_lambda(rhs,lambda_->set_quantity());
932
1167
 
933
1168
    // compute nodal atomic power
934
1169
    compute_nodal_lambda_force(dt);
943
1178
  {
944
1179
    // apply time filtered lambda force
945
1180
    DENS_MAT lambdaZero(nNodes_,nsd_);
946
 
    timeFilter_->apply_pre_step1(lambdaForceFiltered_,lambdaZero,dt);
947
 
    timeFilter_->apply_pre_step1(nodalGhostForceFiltered_,nodalGhostForce_,dt);
 
1181
    timeFilter_->apply_pre_step1(lambdaForceFiltered_->set_quantity(),lambdaZero,dt);
 
1182
    if (nodalGhostForce_) {
 
1183
      timeFilter_->apply_pre_step1(nodalGhostForceFiltered_.set_quantity(),
 
1184
                                   nodalGhostForce_->quantity(),dt);
 
1185
    }
948
1186
  }
949
1187
 
950
1188
  //--------------------------------------------------------
958
1196
    // form rhs :  \int N_I r dV - \sum_g N_Ig^* f_g
959
1197
    // sources are set in ATC transfer
960
1198
    rhs.reset(nNodes_,nsd_);
961
 
    rhs = momentumSource_ - nodalGhostForce_;
962
 
 
 
1199
    rhs = momentumSource_.quantity();
 
1200
    if (nodalGhostForce_) {
 
1201
      rhs -= nodalGhostForce_->quantity();
 
1202
    }
 
1203
    
963
1204
    // (b) for ess. bcs
964
 
    // NOTE can add check to see if set is empty for faster performance
 
1205
    
965
1206
    // form rhs : {sum_a (N_Ia * f_ia) - M_md * (ddupsilon/dt)_I}
966
 
    DENS_MAT rhsPrescribed = -1.*nodalForce_;
967
 
    atcTransfer_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
968
 
    rhsPrescribed = mdMassMatrix_*rhsPrescribed;
969
 
    rhsPrescribed += nodalAtomicForce_;
 
1207
    DENS_MAT rhsPrescribed = -1.*nodalForce_.quantity();
 
1208
    atc_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
 
1209
    rhsPrescribed = (mdMassMatrix_.quantity())*rhsPrescribed;
 
1210
    rhsPrescribed += nodalAtomicForce_->quantity();
970
1211
 
971
1212
    set<pair<int,int> >::const_iterator iter;
972
1213
    for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
975
1216
  }
976
1217
 
977
1218
  //--------------------------------------------------------
978
 
  //  compute_lambda_force
979
 
  //            computes the force induced by lambda
980
 
  //            on the atoms
981
 
  //--------------------------------------------------------
982
 
  void StressFlux::compute_lambda_force(DENS_MAT & lambdaForce)
983
 
  {
984
 
    if (nLocal_>0) {
985
 
      // prolongation to (unique) nodes
986
 
      lambdaForce.reset(nLocal_,nsd_);
987
 
      atcTransfer_->prolong(lambda_,lambdaForce);
988
 
      lambdaForce *= -1.;
989
 
    }
990
 
  }
991
 
 
992
 
  //--------------------------------------------------------
993
1219
  //  compute_nodal_lambda_force
994
1220
  //            computes the force induced on the FE
995
1221
  //            by applying lambdaForce on the atoms
996
1222
  //--------------------------------------------------------
997
1223
  void StressFlux::compute_nodal_lambda_force(double dt)
998
1224
  {
999
 
    atcTransfer_->restrict_volumetric_quantity(lambdaForce_,nodalAtomicLambdaForce_);
 
1225
    DENS_MAT myNodalAtomicLambdaForce = nodalAtomicLambdaForce_->quantity();
1000
1226
    set<pair<int,int> >::const_iterator iter;
1001
1227
    for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
1002
 
      nodalAtomicLambdaForce_(iter->first,iter->second) = 0.;
 
1228
      myNodalAtomicLambdaForce(iter->first,iter->second) = 0.;
1003
1229
    }
1004
1230
 
1005
 
    timeFilter_->apply_post_step1(lambdaForceFiltered_,nodalAtomicLambdaForce_,dt);
 
1231
    timeFilter_->apply_post_step1(lambdaForceFiltered_->set_quantity(),
 
1232
                                  myNodalAtomicLambdaForce,dt);
1006
1233
  }
1007
1234
 
1008
1235
  //--------------------------------------------------------
1013
1240
  //--------------------------------------------------------
1014
1241
  void StressFlux::add_to_rhs(FIELDS & rhs)
1015
1242
  {
1016
 
    rhs[VELOCITY] += nodalAtomicLambdaForce_ + boundaryFlux_[VELOCITY];
1017
 
  }
 
1243
    // compute the kinetostat force
 
1244
    compute_kinetostat(atc_->dt());
1018
1245
 
1019
 
  //--------------------------------------------------------
1020
 
  //  compute_ghost_force
1021
 
  //            computes computes the restricted force on
1022
 
  //            the ghost atoms
1023
 
  //--------------------------------------------------------
1024
 
  void StressFlux::compute_ghost_force(DENS_MAT & nodalGhostForce)
1025
 
  { 
1026
 
    DENS_MAT nodalGhostForceLocal(nNodes_,nsd_);
1027
 
    nodalGhostForce.reset(nNodes_,nsd_);
1028
 
    
1029
 
    if (nLocalGhost_ > 0) {
1030
 
      DENS_MAT ghostForce(nLocalGhost_,nsd_);
1031
 
      for (int i = 0; i < nLocalGhost_; i++) {
1032
 
        int atomIndex = ghostToAtom_(i);
1033
 
        for (int j = 0; j < nsd_; j++) {
1034
 
          ghostForce(i,j) = f_[atomIndex][j];
1035
 
#if false
1036
 
          ghostForce(i,j) -= f0_[atomIndex][j];
1037
 
#endif
1038
 
        }
1039
 
      }
1040
 
      nodalGhostForceLocal = shapeFunctionGhost_.transMat(ghostForce);
1041
 
    }
1042
 
    
1043
 
    LammpsInterface::instance()->allsum(nodalGhostForceLocal.get_ptr(),
1044
 
                                        nodalGhostForce.get_ptr(), nodalGhostForce.size());
1045
 
    // convert from Lammps force units to ATC force units
1046
 
    nodalGhostForce *= LammpsInterface::instance()->ftm2v();
 
1246
    rhs[VELOCITY] += nodalAtomicLambdaForce_->quantity() + boundaryFlux_[VELOCITY].quantity();
1047
1247
  }
1048
1248
 
1049
1249
  //--------------------------------------------------------
1051
1251
  //            uses existing lambda to modify given
1052
1252
  //            atomic quantity
1053
1253
  //--------------------------------------------------------
1054
 
  void StressFlux::apply_to_atoms(double ** atomicVelocity,
 
1254
  void StressFlux::apply_to_atoms(PerAtomQuantity<double> * atomVelocities,
1055
1255
                                  const DENS_MAT & lambdaForce,
1056
1256
                                  double dt)
1057
1257
  {
1058
 
    if (nLocal_>0) {
1059
 
      // explicit update
1060
 
      for (int i = 0; i < nLocal_; i++) {
1061
 
        for (int j = 0; j < nsd_; j++) {
1062
 
          atomicVelocity[internalToAtom_(i)][j] += dt*lambdaForce(i,j)/atomicMass_(i);
1063
 
        }
1064
 
      }
1065
 
    }
 
1258
    _deltaVelocity_ = lambdaForce;
 
1259
    _deltaVelocity_ /= atomMasses_->quantity();
 
1260
    _deltaVelocity_ *= dt;
 
1261
    *atomVelocities += _deltaVelocity_;
1066
1262
  }
1067
1263
 
1068
1264
  //--------------------------------------------------------
1072
1268
  //--------------------------------------------------------
1073
1269
  void StressFlux::reset_filtered_ghost_force(DENS_MAT & target)
1074
1270
  {
1075
 
    for (int i = 0; i < nNodes_; ++i)
1076
 
      for (int j = 0; j < nsd_; ++j)
1077
 
        nodalGhostForceFiltered_(i,j) = target(i,j);
 
1271
    nodalGhostForceFiltered_ = target;
1078
1272
  }
1079
1273
 
1080
1274
  //--------------------------------------------------------
1081
1275
  //  output:
1082
1276
  //    adds all relevant output to outputData
1083
1277
  //--------------------------------------------------------
1084
 
  void StressFlux::output(double dt, OUTPUT_LIST & outputData)
1085
 
  {
1086
 
    outputData["lambda"] = &(lambda_);
1087
 
    outputData["nodalLambdaForce"] = &(nodalAtomicLambdaForce_);
 
1278
  void StressFlux::output(OUTPUT_LIST & outputData)
 
1279
  {
 
1280
    _nodalAtomicLambdaForceOut_ = nodalAtomicLambdaForce_->quantity();
 
1281
    DENS_MAT & lambda(lambda_->set_quantity());
 
1282
    if ((atc_->lammps_interface())->rank_zero()) {
 
1283
      outputData[regulatorPrefix_+"Lambda"] = &lambda;
 
1284
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_nodalAtomicLambdaForceOut_);
 
1285
    }
 
1286
  }
 
1287
 
 
1288
  //--------------------------------------------------------
 
1289
  //--------------------------------------------------------
 
1290
  //  Class StressFluxGhost
 
1291
  //--------------------------------------------------------
 
1292
  //--------------------------------------------------------
 
1293
 
 
1294
  //--------------------------------------------------------
 
1295
  //  Constructor
 
1296
  //         Grab references to ATC and kinetostat data
 
1297
  //--------------------------------------------------------
 
1298
  StressFluxGhost::StressFluxGhost(Kinetostat * kinetostat) :
 
1299
    StressFlux(kinetostat)
 
1300
  {
 
1301
    // flag for performing boundary flux calculation
 
1302
    fieldMask_(VELOCITY,FLUX) = false;
 
1303
  }
 
1304
 
 
1305
  //--------------------------------------------------------
 
1306
  //  constructor_transfers
 
1307
  //    instantiates or obtains all dependency managed data
 
1308
  //--------------------------------------------------------
 
1309
  void StressFluxGhost::construct_transfers()
 
1310
  {
 
1311
    StressFlux::construct_transfers();
 
1312
    if (!nodalGhostForce_) {
 
1313
      throw ATC_Error("StressFluxGhost::StressFluxGhost - ghost atoms must be specified");
 
1314
    }
 
1315
  }
 
1316
 
 
1317
  //--------------------------------------------------------
 
1318
  //  compute_boundary_flux:
 
1319
  //    computes the boundary flux to be consistent with
 
1320
  //    the controller
 
1321
  //--------------------------------------------------------
 
1322
  void StressFluxGhost::compute_boundary_flux(FIELDS & fields)
 
1323
  {
 
1324
    // This is only used in computation of atomic sources
 
1325
    boundaryFlux_[VELOCITY] = 0.;
 
1326
  }
 
1327
 
 
1328
  //--------------------------------------------------------
 
1329
  //  add_to_rhs:
 
1330
  //            determines what if any contributions to the
 
1331
  //            finite element equations are needed for
 
1332
  //            consistency with the kinetostat
 
1333
  //--------------------------------------------------------
 
1334
  void StressFluxGhost::add_to_rhs(FIELDS & rhs)
 
1335
  {
 
1336
    // compute the kinetostat force
 
1337
    compute_kinetostat(atc_->dt());
 
1338
 
 
1339
    // uses ghost force as the boundary flux to add to the RHS
 
1340
    rhs[VELOCITY] += nodalAtomicLambdaForce_->quantity() + nodalGhostForce_->quantity();
 
1341
  }
 
1342
 
 
1343
    //--------------------------------------------------------
 
1344
  //  set_kinetostat_rhs
 
1345
  //            sets up the RHS of the kinetostat equations
 
1346
  //            for the coupling parameter lambda
 
1347
  //--------------------------------------------------------
 
1348
  void StressFluxGhost::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
 
1349
  {
 
1350
    // (a) for flux based : 
 
1351
    // form rhs :  \int N_I r dV - \sum_g N_Ig^* f_g
 
1352
    // sources are set in ATC transfer
 
1353
    rhs.reset(nNodes_,nsd_);
 
1354
    rhs = momentumSource_.quantity();
 
1355
    
 
1356
    // (b) for ess. bcs
 
1357
    
 
1358
    // form rhs : {sum_a (N_Ia * f_ia) - M_md * (ddupsilon/dt)_I}
 
1359
    DENS_MAT rhsPrescribed = -1.*nodalForce_.quantity();
 
1360
    atc_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
 
1361
    rhsPrescribed = (mdMassMatrix_.quantity())*rhsPrescribed;
 
1362
    rhsPrescribed += nodalAtomicForce_->quantity();
 
1363
 
 
1364
    set<pair<int,int> >::const_iterator iter;
 
1365
    for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
 
1366
      rhs(iter->first,iter->second) = rhsPrescribed(iter->first,iter->second);
 
1367
    }
1088
1368
  }
1089
1369
 
1090
1370
  //--------------------------------------------------------
1098
1378
  //--------------------------------------------------------
1099
1379
  StressFluxFiltered::StressFluxFiltered(Kinetostat * kinetostat) :
1100
1380
    StressFlux(kinetostat),
1101
 
    nodalAtomicVelocity_(atcTransfer_->get_atomic_field(VELOCITY))
 
1381
    nodalAtomicVelocity_(atc_->nodal_atomic_field(VELOCITY))
1102
1382
  {
1103
1383
    // do nothing
1104
1384
  }
1105
 
 
1106
 
#if false
1107
 
  //--------------------------------------------------------
1108
 
  //  apply:
1109
 
  //    apply the kinetostat to the atoms
1110
 
  //--------------------------------------------------------
1111
 
  void StressFluxFiltered::apply_pre_corrector(double dt)
1112
 
  {
1113
 
    // NOTE currently not implemented, below based on GlcVelocityFiltered
1114
 
 
1115
 
    // get neccesary data
1116
 
    double coef = 1./(timeFilter_->get_unfiltered_coefficient_pre_s1(dt));
1117
 
   
1118
 
    // form rhs : sum_a (hatN_Ia * x_ai) - (\dot{Upsilon})_Ii
1119
 
    DENS_MAT rhs(nNodes_,nsd_);
1120
 
    for (int i = 0; i < nNodes_; i++)
1121
 
      for (int j = 0; j < nsd_; j++)
1122
 
        rhs(i,j) += coef*(nodalAtomicVelocities_(i,j) - nodalVelocities_(i,j));
1123
 
    DENS_MAT rhsOverlap(nNodeOverlap_,nsd_);
1124
 
    atcTransfer_->map_unique_to_overlap(rhs, rhsOverlap);
1125
 
    
1126
 
    // solve matrix equation and map back to all nodes
1127
 
    DENS_MAT lambdaOverlap(nNodeOverlap_,nsd_);
1128
 
    DIAG_MAT weights;
1129
 
    set_weights(weights);
1130
 
    solve_for_lambda(rhsOverlap, weights, lambdaOverlap);
1131
 
    atcTransfer_->map_overlap_to_unique(lambdaOverlap,lambda_);
1132
 
        
1133
 
    // apply lambda to the atoms and nodal atomic fields
1134
 
    DENS_MAT lambdaRestricted(nNodes_,nsd_);
1135
 
    apply_lambda_to_atoms(v_,lambdaRestricted);
1136
 
    timeFilter_->apply_post_step1(nodalAtomicVelocities_,-1.*lambdaRestricted,dt);
1137
 
    timeFilter_->apply_post_step1(nodalLambdaForce_,(-1./dt)*lambdaRestricted,dt);
1138
 
    DENS_MAT nodalLambdaRoc(nodalLambdaForce_.nRows,nodalLambdaForce_.nCols());
1139
 
    atcTransfer_->apply_inverse_mass_matrix(nodalLambdaForce_,
1140
 
                                            nodalLambdaRoc,
1141
 
                                            VELOCITY);
1142
 
    atcTransfer_->apply_internal_atomic_source(nodalVelocities_,
1143
 
                                               nodalLambdaRoc,
1144
 
                                               dt);
1145
 
 
1146
 
  }
1147
 
#endif
1148
 
//   //--------------------------------------------------------
1149
 
//   //  apply_pre_filtering
1150
 
//   //            applies first step of filtering to
1151
 
//   //            relevant variables
1152
 
//   //--------------------------------------------------------
1153
 
//   void StressFluxFiltered::apply_pre_filtering(double dt)
1154
 
//   {
1155
 
//     // apply time filtered lambda to atomic fields
1156
 
//     StressFlux::apply_pre_filtering(dt);
1157
 
//     nodalAtomicForce_ += lambdaForceFiltered;
1158
 
//   }
1159
1385
 
1160
1386
  //--------------------------------------------------------
1161
1387
  //  set_kinetostat_rhs
1169
1395
    // form rhs :  \int N_I r dV - \sum_g N_Ig^* f_g
1170
1396
    // sources are set in ATC transfer
1171
1397
    rhs.reset(nNodes_,nsd_);
1172
 
    rhs = momentumSource_ - nodalGhostForceFiltered_;
 
1398
    rhs = momentumSource_.quantity() - nodalGhostForceFiltered_.quantity();
1173
1399
    
1174
1400
    // (b) for ess. bcs
1175
 
    // NOTE can add check to see if set is empty for faster performance
 
1401
    
1176
1402
    // form rhs : {sum_a (N_Ia * f_ia) - M_md * (ddupsilon/dt)_I}
1177
 
    DENS_MAT rhsPrescribed = -1.*nodalForce_;
1178
 
    atcTransfer_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
1179
 
    rhsPrescribed = mdMassMatrix_*rhsPrescribed;
1180
 
    rhsPrescribed += nodalAtomicForce_;
 
1403
    DENS_MAT rhsPrescribed = -1.*nodalForce_.quantity();
 
1404
    atc_->apply_inverse_mass_matrix(rhsPrescribed,VELOCITY);
 
1405
    rhsPrescribed = (mdMassMatrix_.quantity())*rhsPrescribed;
 
1406
    rhsPrescribed += nodalAtomicForce_->quantity();
1181
1407
 
1182
1408
    set<pair<int,int> >::const_iterator iter;
1183
1409
    for (iter = hooverNodes_.begin(); iter != hooverNodes_.end(); ++iter) {
1185
1411
    }
1186
1412
 
1187
1413
    // adjust for application of current lambda force
1188
 
    rhs += lambdaForceFiltered_;
 
1414
    rhs += lambdaForceFiltered_->quantity();
1189
1415
 
1190
1416
    // correct for time filtering
1191
 
    rhs *= 1./(timeFilter_->get_unfiltered_coefficient_pre_s1(dt));
 
1417
    rhs *= 1./(timeFilter_->unfiltered_coefficient_pre_s1(dt));
1192
1418
  }
1193
1419
 
1194
1420
  //--------------------------------------------------------
1196
1422
  //            uses existing lambda to modify given
1197
1423
  //            atomic quantity
1198
1424
  //--------------------------------------------------------
1199
 
  void StressFluxFiltered::apply_to_atoms(double ** atomicVelocity,
 
1425
  void StressFluxFiltered::apply_to_atoms(PerAtomQuantity<double> * atomVelocities,
1200
1426
                                          const DENS_MAT & lambdaForce,
1201
1427
                                          double dt)
1202
1428
  {
1203
 
    StressFlux::apply_to_atoms(atomicVelocity,lambdaForce,dt);
 
1429
    StressFlux::apply_to_atoms(atomVelocities,lambdaForce,dt);
1204
1430
 
1205
1431
    // add in corrections to filtered nodal atomice velocity
1206
1432
    DENS_MAT velocityRoc(nNodes_,nsd_);
1207
 
    atcTransfer_->apply_inverse_md_mass_matrix(lambdaForceFiltered_,
1208
 
                                               velocityRoc,
1209
 
                                               VELOCITY);
 
1433
    atc_->apply_inverse_md_mass_matrix(lambdaForceFiltered_->quantity(),
 
1434
                                       velocityRoc,
 
1435
                                       VELOCITY);
1210
1436
    nodalAtomicVelocity_ += dt*velocityRoc;
1211
1437
  }
1212
1438
 
1218
1444
  //--------------------------------------------------------
1219
1445
  void StressFluxFiltered::add_to_rhs(FIELDS & rhs)
1220
1446
  {
1221
 
    rhs[VELOCITY] += lambdaForceFiltered_ + boundaryFlux_[VELOCITY];
1222
 
  }
1223
 
 
1224
 
  //--------------------------------------------------------
1225
 
  //  output:
1226
 
  //    adds all relevant output to outputData
1227
 
  //--------------------------------------------------------
1228
 
  void StressFluxFiltered::output(double dt, OUTPUT_LIST & outputData)
1229
 
  {
1230
 
    outputData["lambda"] = &(lambda_);
1231
 
    outputData["nodalLambdaForce"] = &(lambdaForceFiltered_);
1232
 
  }
1233
 
  
 
1447
    // compute kinetostat forces
 
1448
    compute_kinetostat(atc_->dt());
 
1449
    rhs[VELOCITY] += lambdaForceFiltered_->quantity() + boundaryFlux_[VELOCITY].quantity();
 
1450
  }
 
1451
 
 
1452
  //--------------------------------------------------------
 
1453
  //  output:
 
1454
  //    adds all relevant output to outputData
 
1455
  //--------------------------------------------------------
 
1456
  void StressFluxFiltered::output(OUTPUT_LIST & outputData)
 
1457
  {
 
1458
    DENS_MAT & lambda(lambda_->set_quantity());
 
1459
    DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
 
1460
    if ((atc_->lammps_interface())->rank_zero()) {
 
1461
      outputData[regulatorPrefix_+"Lambda"] = &lambda;
 
1462
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &lambdaForceFiltered;
 
1463
    }
 
1464
  }
 
1465
 
 
1466
  //--------------------------------------------------------
 
1467
  //--------------------------------------------------------
 
1468
  //  Class KinetostatGlcFs
 
1469
  //--------------------------------------------------------
 
1470
  //--------------------------------------------------------
 
1471
 
 
1472
  //--------------------------------------------------------
 
1473
  //  Constructor
 
1474
  //         Grab references to ATC and kinetostat data
 
1475
  //--------------------------------------------------------
 
1476
  KinetostatGlcFs::KinetostatGlcFs(Kinetostat * kinetostat,
 
1477
                                   const string & regulatorPrefix) :
 
1478
    KinetostatShapeFunction(kinetostat,regulatorPrefix),
 
1479
    velocity_(atc_->field(VELOCITY))
 
1480
  {
 
1481
    // constuct/obtain data corresponding to stage 3 of ATC_Method::initialize
 
1482
    nodalAtomicLambdaForce_ = kinetostat_->regulator_data(regulatorPrefix_+"NodalAtomicLambdaForce",nsd_);
 
1483
  }
 
1484
 
 
1485
  //--------------------------------------------------------
 
1486
  //  constructor_transfers
 
1487
  //    instantiates or obtains all dependency managed data
 
1488
  //--------------------------------------------------------
 
1489
  void KinetostatGlcFs::construct_transfers()
 
1490
  {
 
1491
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
1492
 
 
1493
    // base class transfers
 
1494
    KinetostatShapeFunction::construct_transfers();
 
1495
 
 
1496
    // get data from manager
 
1497
    nodalAtomicMomentum_ = interscaleManager.dense_matrix("NodalAtomicMomentum");
 
1498
 
 
1499
    // atomic force induced by kinetostat
 
1500
    PerAtomQuantity<double> * atomLambdas = interscaleManager.per_atom_quantity(regulatorPrefix_+"AtomLambdaMomentum");
 
1501
    atomKinetostatForce_ = new AtomicKinetostatForceStress(atc_,atomLambdas);
 
1502
    interscaleManager.add_per_atom_quantity(atomKinetostatForce_,
 
1503
                                            regulatorPrefix_+"AtomKinetostatForce");
 
1504
  }
 
1505
 
 
1506
  //--------------------------------------------------------
 
1507
  //  initialize
 
1508
  //    initializes all method data
 
1509
  //--------------------------------------------------------
 
1510
  void KinetostatGlcFs::initialize()
 
1511
  {
 
1512
    KinetostatShapeFunction::initialize();
 
1513
 
 
1514
    // set up workspaces
 
1515
    _deltaMomentum_.reset(nNodes_,nsd_);
 
1516
    _lambdaForceOutput_.reset(nNodes_,nsd_);
 
1517
    
 
1518
    TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
 
1519
    if (!timeFilterManager->end_equilibrate()) {
 
1520
      // we should reset lambda and lambdaForce to zero in this case
 
1521
      // implies an initial condition of 0 for the filtered nodal lambda power
 
1522
      // initial conditions will always be needed when using time filtering
 
1523
      // however, the fractional step scheme must assume the instantaneous
 
1524
      // nodal lambda power is 0 initially because all quantities are in delta form
 
1525
      *lambda_ = 0.; // ensures initial lambda force is zero
 
1526
      *nodalAtomicLambdaForce_ = 0.; // momentum change due to kinetostat
 
1527
      *lambdaForceFiltered_ = 0.; // filtered momentum change due to kinetostats
 
1528
    }
 
1529
    else {
 
1530
      // we can grab lambda power variables using time integrator and atc transfer in cases for equilibration
 
1531
    }
 
1532
 
 
1533
    // sets up time filter for cases where variables temporally filtered
 
1534
    if (timeFilterManager->need_reset()) {
 
1535
      // the form of this integrator implies no time filters that require history data can be used
 
1536
      timeFilter_->initialize(nodalAtomicLambdaForce_->quantity());
 
1537
    }
 
1538
 
 
1539
    compute_rhs_map();
 
1540
  }
 
1541
 
 
1542
  //--------------------------------------------------------
 
1543
  //  compute_rhs_map
 
1544
  //    creates mapping from all nodes to those to which
 
1545
  //    the kinetostat applies
 
1546
  //--------------------------------------------------------
 
1547
  void KinetostatGlcFs::compute_rhs_map()
 
1548
  {
 
1549
    rhsMap_.resize(overlapToNodeMap_->nRows(),1);
 
1550
    DENS_MAT rhsMapGlobal(nNodes_,1);
 
1551
    const set<int> & applicationNodes(applicationNodes_->quantity());
 
1552
    for (int i = 0; i < nNodes_; i++) {
 
1553
      if (applicationNodes.find(i) != applicationNodes.end()) {
 
1554
        rhsMapGlobal(i,0) = 1.;
 
1555
      }
 
1556
      else {
 
1557
        rhsMapGlobal(i,0) = 0.;
 
1558
      }
 
1559
    }
 
1560
    map_unique_to_overlap(rhsMapGlobal,rhsMap_);
 
1561
  }
 
1562
 
 
1563
  //--------------------------------------------------------
 
1564
  //  apply_pre_predictor:
 
1565
  //    apply the kinetostat to the atoms in the
 
1566
  //    pre-predictor integration phase
 
1567
  //--------------------------------------------------------
 
1568
  void KinetostatGlcFs::apply_pre_predictor(double dt)
 
1569
  {
 
1570
    DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
 
1571
    DENS_MAT & nodalAtomicLambdaForce(nodalAtomicLambdaForce_->set_quantity());
 
1572
 
 
1573
    // update filtered forces
 
1574
    timeFilter_->apply_pre_step1(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
 
1575
 
 
1576
    // apply lambda force to atoms and compute instantaneous lambda force 
 
1577
    this->apply_to_atoms(atomVelocities_,nodalAtomicMomentum_,
 
1578
                         atomKinetostatForce_->quantity(),
 
1579
                         nodalAtomicLambdaForce,0.5*dt);
 
1580
 
 
1581
    // update nodal variables for first half of timestep
 
1582
    this->add_to_momentum(nodalAtomicLambdaForce,_deltaMomentum_,0.5*dt);
 
1583
    atc_->apply_inverse_mass_matrix(_deltaMomentum_,VELOCITY);
 
1584
    velocity_ += _deltaMomentum_;
 
1585
 
 
1586
    // start update of filtered lambda force
 
1587
    nodalAtomicLambdaForce = 0.;
 
1588
    timeFilter_->apply_post_step1(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
 
1589
  }
 
1590
 
 
1591
  //--------------------------------------------------------
 
1592
  //  apply_post_corrector:
 
1593
  //    apply the kinetostat to the atoms in the
 
1594
  //    post-corrector integration phase
 
1595
  //--------------------------------------------------------
 
1596
  void KinetostatGlcFs::apply_post_corrector(double dt)
 
1597
  {
 
1598
    // compute the kinetostat equation and update lambda
 
1599
    this->compute_lambda(dt);
 
1600
 
 
1601
    DENS_MAT & lambdaForceFiltered(lambdaForceFiltered_->set_quantity());
 
1602
    DENS_MAT & nodalAtomicLambdaForce(nodalAtomicLambdaForce_->set_quantity());
 
1603
 
 
1604
    // update filtered force
 
1605
    timeFilter_->apply_pre_step1(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
 
1606
 
 
1607
    // apply lambda force to atoms and compute instantaneous lambda force 
 
1608
    this->apply_to_atoms(atomVelocities_,nodalAtomicMomentum_,
 
1609
                         atomKinetostatForce_->quantity(),
 
1610
                         nodalAtomicLambdaForce,0.5*dt);
 
1611
 
 
1612
    // update nodal variables for first half of timestep
 
1613
    this->add_to_momentum(nodalAtomicLambdaForce,_deltaMomentum_,0.5*dt);
 
1614
    nodalAtomicLambdaForce *= 2./dt;
 
1615
    atc_->apply_inverse_mass_matrix(_deltaMomentum_,VELOCITY);
 
1616
    velocity_ += _deltaMomentum_;
 
1617
 
 
1618
    // start update of filtered lambda force
 
1619
    timeFilter_->apply_post_step2(lambdaForceFiltered,nodalAtomicLambdaForce,dt);
 
1620
  }
 
1621
 
 
1622
  //--------------------------------------------------------
 
1623
  //  compute_kinetostat
 
1624
  //            manages the solution and application of the
 
1625
  //            kinetostat equations and variables
 
1626
  //--------------------------------------------------------
 
1627
  void KinetostatGlcFs::compute_lambda(double dt)
 
1628
  {
 
1629
    // set up rhs for lambda equation
 
1630
    this->set_kinetostat_rhs(rhs_,0.5*dt);
 
1631
 
 
1632
    // solve linear system for lambda
 
1633
    DENS_MAT & lambda(lambda_->set_quantity());
 
1634
    solve_for_lambda(rhs_,lambda);
 
1635
  }
 
1636
 
 
1637
  //--------------------------------------------------------
 
1638
  //  apply_lambda_to_atoms
 
1639
  //            uses existing lambda to modify given
 
1640
  //            atomic quantity
 
1641
  //--------------------------------------------------------
 
1642
  void KinetostatGlcFs::apply_to_atoms(PerAtomQuantity<double> * atomVelocity,
 
1643
                                      const DENS_MAN * nodalAtomicMomentum,
 
1644
                                      const DENS_MAT & lambdaForce,
 
1645
                                      DENS_MAT & nodalAtomicLambdaForce,
 
1646
                                      double dt)
 
1647
  {
 
1648
    // compute initial contributions to lambda force
 
1649
    nodalAtomicLambdaForce = nodalAtomicMomentum->quantity();
 
1650
    nodalAtomicLambdaForce *= -1.;
 
1651
 
 
1652
    // apply lambda force to atoms
 
1653
    _velocityDelta_ = lambdaForce;
 
1654
    _velocityDelta_ /= atomMasses_->quantity();
 
1655
    _velocityDelta_ *= dt;
 
1656
    (*atomVelocity) += _velocityDelta_;
 
1657
 
 
1658
    // finalize lambda force
 
1659
    nodalAtomicLambdaForce += nodalAtomicMomentum->quantity();
 
1660
  }
 
1661
 
 
1662
  //--------------------------------------------------------
 
1663
  //  output:
 
1664
  //    adds all relevant output to outputData
 
1665
  //--------------------------------------------------------
 
1666
  void KinetostatGlcFs::output(OUTPUT_LIST & outputData)
 
1667
  {
 
1668
    _lambdaForceOutput_ = nodalAtomicLambdaForce_->quantity();
 
1669
    // approximate value for lambda force
 
1670
    double dt =  LammpsInterface::instance()->dt();
 
1671
    _lambdaForceOutput_ *= (2./dt);
 
1672
    DENS_MAT & lambda(lambda_->set_quantity());
 
1673
    if ((atc_->lammps_interface())->rank_zero()) {
 
1674
      outputData[regulatorPrefix_+"Lambda"] = &lambda;
 
1675
      outputData[regulatorPrefix_+"NodalLambdaForce"] = &(_lambdaForceOutput_);
 
1676
    }
 
1677
  }
 
1678
 
 
1679
  //--------------------------------------------------------
 
1680
  //--------------------------------------------------------
 
1681
  //  Class KinetostatFlux
 
1682
  //--------------------------------------------------------
 
1683
  //--------------------------------------------------------
 
1684
 
 
1685
  //--------------------------------------------------------
 
1686
  //  Constructor
 
1687
  //         Grab references to ATC and kinetostat data
 
1688
  //--------------------------------------------------------
 
1689
  KinetostatFlux::KinetostatFlux(Kinetostat * kinetostat,
 
1690
                                 const string & regulatorPrefix) :
 
1691
    KinetostatGlcFs(kinetostat,regulatorPrefix),
 
1692
    momentumSource_(atc_->atomic_source(VELOCITY)),
 
1693
    nodalGhostForce_(NULL),
 
1694
    nodalGhostForceFiltered_(NULL)
 
1695
  {
 
1696
    // flag for performing boundary flux calculation
 
1697
    fieldMask_(VELOCITY,FLUX) = true;
 
1698
 
 
1699
    // constuct/obtain data corresponding to stage 3 of ATC_Method::initialize
 
1700
    nodalGhostForceFiltered_ = kinetostat_->regulator_data(regulatorPrefix_+"NodalGhostForceFiltered",nsd_);
 
1701
  }
 
1702
 
 
1703
  //--------------------------------------------------------
 
1704
  //  constructor_transfers
 
1705
  //    instantiates or obtains all dependency managed data
 
1706
  //--------------------------------------------------------
 
1707
  void KinetostatFlux::construct_transfers()
 
1708
  {
 
1709
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
1710
 
 
1711
    // set up node mappings
 
1712
    create_node_maps();
 
1713
 
 
1714
    
 
1715
    // set up linear solver
 
1716
    // set up data for linear solver
 
1717
    shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
 
1718
    interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
 
1719
                                                 regulatorPrefix_+"LambdaCouplingMatrixMomentum");
 
1720
    if (elementMask_) {
 
1721
      lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
 
1722
      interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
 
1723
                                                  regulatorPrefix_+"LambdaAtomMap");
 
1724
    }
 
1725
    if (atomicRegulator_->use_localized_lambda()) {
 
1726
      linearSolverType_ = AtomicRegulator::RSL_SOLVE;
 
1727
    }
 
1728
    else {
 
1729
      linearSolverType_ = AtomicRegulator::CG_SOLVE;
 
1730
    }
 
1731
 
 
1732
    // base class transfers
 
1733
    KinetostatGlcFs::construct_transfers();
 
1734
 
 
1735
    // sets up space for ghost force related variables
 
1736
    if (atc_->groupbit_ghost()) {
 
1737
      MatrixDependencyManager<DenseMatrix, int> * nodeToOverlapMap = 
 
1738
        interscaleManager.dense_matrix_int(regulatorPrefix_+"NodeToOverlapMap");
 
1739
      GhostCouplingMatrix * shapeFunctionGhost = new GhostCouplingMatrix(atc_,interscaleManager.per_atom_sparse_matrix("InterpolantGhost"),
 
1740
                                                                         regulatedNodes_,
 
1741
                                                                         nodeToOverlapMap);
 
1742
      interscaleManager.add_sparse_matrix(shapeFunctionGhost,
 
1743
                                          regulatorPrefix_+"GhostCouplingMatrix");
 
1744
      FundamentalAtomQuantity * atomGhostForce = interscaleManager.fundamental_atom_quantity(LammpsInterface::ATOM_FORCE,
 
1745
                                                                                             GHOST);
 
1746
      nodalGhostForce_ = new AtfShapeFunctionRestriction(atc_,atomGhostForce,
 
1747
                                                         shapeFunctionGhost);
 
1748
      interscaleManager.add_dense_matrix(nodalGhostForce_,
 
1749
                                         regulatorPrefix_+"NodalGhostForce");
 
1750
    }
 
1751
  }
 
1752
 
 
1753
  //--------------------------------------------------------
 
1754
  //  initialize
 
1755
  //    initializes all method data
 
1756
  //--------------------------------------------------------
 
1757
  void KinetostatFlux::initialize()
 
1758
  {
 
1759
    KinetostatGlcFs::initialize();
 
1760
    
 
1761
    TimeFilterManager * timeFilterManager = atc_->time_filter_manager();
 
1762
    if (!timeFilterManager->end_equilibrate()) {
 
1763
      // we should reset lambda and lambdaForce to zero in this case
 
1764
      // implies an initial condition of 0 for the filtered nodal lambda power
 
1765
      // initial conditions will always be needed when using time filtering
 
1766
      // however, the fractional step scheme must assume the instantaneous
 
1767
      // nodal lambda power is 0 initially because all quantities are in delta form
 
1768
      *nodalGhostForceFiltered_ = 0.; // filtered force from ghost atoms
 
1769
    }
 
1770
    else {
 
1771
      // we can grab lambda power variables using time integrator and atc transfer in cases for equilibration
 
1772
    }
 
1773
  }
 
1774
 
 
1775
  //--------------------------------------------------------
 
1776
  //  construct_regulated_nodes:
 
1777
  //    constructs the set of nodes being regulated
 
1778
  //--------------------------------------------------------
 
1779
  void KinetostatFlux::construct_regulated_nodes()
 
1780
  {
 
1781
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
1782
 
 
1783
    // matrix requires all entries even if localized for correct lumping
 
1784
    regulatedNodes_ = new RegulatedNodes(atc_);
 
1785
    interscaleManager.add_set_int(regulatedNodes_,
 
1786
                                  regulatorPrefix_+"KinetostatRegulatedNodes");
 
1787
 
 
1788
    // if localized monitor nodes with applied fluxes
 
1789
    if (atomicRegulator_->use_localized_lambda()) {
 
1790
      if ((kinetostat_->coupling_mode() == Kinetostat::FLUX) && (atomicRegulator_->boundary_integration_type() != NO_QUADRATURE)) {
 
1791
        // include boundary nodes
 
1792
        applicationNodes_ = new FluxBoundaryNodes(atc_);
 
1793
        
 
1794
        boundaryNodes_ = new BoundaryNodes(atc_);
 
1795
        interscaleManager.add_set_int(boundaryNodes_,
 
1796
                                      regulatorPrefix_+"KinetostatBoundaryNodes");
 
1797
      }
 
1798
      else {
 
1799
        // fluxed nodes only
 
1800
        applicationNodes_ = new FluxNodes(atc_);
 
1801
      }
 
1802
      interscaleManager.add_set_int(applicationNodes_,
 
1803
                                    regulatorPrefix_+"KinetostatApplicationNodes");
 
1804
    }
 
1805
    else {
 
1806
      applicationNodes_ = regulatedNodes_;
 
1807
    }
 
1808
 
 
1809
    // special set of boundary elements for boundary flux quadrature  
 
1810
    if ((atomicRegulator_->boundary_integration_type() == FE_INTERPOLATION)
 
1811
        && (atomicRegulator_->use_localized_lambda())) {
 
1812
      elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
 
1813
      interscaleManager.add_dense_matrix_bool(elementMask_,
 
1814
                                                   regulatorPrefix_+"BoundaryElementMask");
 
1815
    }
 
1816
  }
 
1817
 
 
1818
  //--------------------------------------------------------
 
1819
  //  apply_pre_predictor:
 
1820
  //    apply the kinetostat to the atoms in the
 
1821
  //    pre-predictor integration phase
 
1822
  //--------------------------------------------------------
 
1823
  void KinetostatFlux::apply_pre_predictor(double dt)
 
1824
  {
 
1825
    // update filtered forces
 
1826
    if (nodalGhostForce_) {
 
1827
      timeFilter_->apply_pre_step1(nodalGhostForceFiltered_->set_quantity(),
 
1828
                                   nodalGhostForce_->quantity(),dt);
 
1829
    }
 
1830
 
 
1831
    KinetostatGlcFs::apply_pre_predictor(dt);
 
1832
  }
 
1833
 
 
1834
  //--------------------------------------------------------
 
1835
  //  apply_post_corrector:
 
1836
  //    apply the kinetostat to the atoms in the
 
1837
  //    post-corrector integration phase
 
1838
  //--------------------------------------------------------
 
1839
  void KinetostatFlux::apply_post_corrector(double dt)
 
1840
  {
 
1841
    // update filtered ghost force
 
1842
    if (nodalGhostForce_) {
 
1843
      timeFilter_->apply_post_step1(nodalGhostForceFiltered_->set_quantity(),
 
1844
                                    nodalGhostForce_->quantity(),dt);
 
1845
    }
 
1846
    
 
1847
    // compute the kinetostat equation and update lambda
 
1848
    KinetostatGlcFs::apply_post_corrector(dt);
 
1849
  }
 
1850
 
 
1851
  //--------------------------------------------------------
 
1852
  //  set_kinetostat_rhs
 
1853
  //            sets up the RHS of the kinetostat equations
 
1854
  //            for the coupling parameter lambda
 
1855
  //--------------------------------------------------------
 
1856
  void KinetostatFlux::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
 
1857
  {
 
1858
    // (a) for flux based : 
 
1859
    // form rhs :  \int N_I r dV - \sum_g N_Ig^* f_g
 
1860
    // sources are set in ATC transfer
 
1861
    rhs.reset(nNodes_,nsd_);
 
1862
    const DENS_MAT & momentumSource(momentumSource_.quantity());
 
1863
    const set<int> & applicationNodes(applicationNodes_->quantity());
 
1864
    set<int>::const_iterator iNode;
 
1865
    for (iNode = applicationNodes.begin(); iNode != applicationNodes.end(); iNode++) {
 
1866
      for (int j = 0; j < nsd_; j++) {
 
1867
        rhs(*iNode,j) = momentumSource(*iNode,j);
 
1868
      }
 
1869
    }
 
1870
 
 
1871
    // add ghost forces, if needed
 
1872
    if (nodalGhostForce_) {
 
1873
      const DENS_MAT & nodalGhostForce(nodalGhostForce_->quantity());
 
1874
      for (iNode = applicationNodes.begin(); iNode != applicationNodes.end(); iNode++) {
 
1875
        for (int j = 0; j < nsd_; j++) {
 
1876
          rhs(*iNode,j) -= nodalGhostForce(*iNode,j);
 
1877
        }
 
1878
      }
 
1879
    }
 
1880
  }
 
1881
 
 
1882
  //--------------------------------------------------------
 
1883
  //  add_to_momentum:
 
1884
  //            determines what if any contributions to the
 
1885
  //            finite element equations are needed for
 
1886
  //            consistency with the kinetostat
 
1887
  //--------------------------------------------------------
 
1888
  void KinetostatFlux::add_to_momentum(const DENS_MAT & nodalLambdaForce,
 
1889
                                       DENS_MAT & deltaMomentum,
 
1890
                                       double dt)
 
1891
  {
 
1892
    deltaMomentum.resize(nNodes_,nsd_);
 
1893
    const DENS_MAT & boundaryFlux(boundaryFlux_[VELOCITY].quantity());
 
1894
    for (int i = 0; i < nNodes_; i++) {
 
1895
      for (int j = 0; j < nsd_; j++) {
 
1896
        deltaMomentum(i,j) = nodalLambdaForce(i,j) + dt*boundaryFlux(i,j);
 
1897
      }
 
1898
    }
 
1899
  }
 
1900
 
 
1901
  //--------------------------------------------------------
 
1902
  //  reset_filtered_ghost_force:
 
1903
  //    resets the kinetostat generated ghost force to a
 
1904
  //    prescribed value
 
1905
  //--------------------------------------------------------
 
1906
  void KinetostatFlux::reset_filtered_ghost_force(DENS_MAT & target)
 
1907
  {
 
1908
    (*nodalGhostForceFiltered_) = target;
 
1909
  }
 
1910
 
 
1911
  //--------------------------------------------------------
 
1912
  //--------------------------------------------------------
 
1913
  //  Class KinetostatFluxGhost
 
1914
  //--------------------------------------------------------
 
1915
  //--------------------------------------------------------
 
1916
 
 
1917
  //--------------------------------------------------------
 
1918
  //  Constructor
 
1919
  //         Grab references to ATC and kinetostat data
 
1920
  //--------------------------------------------------------
 
1921
  KinetostatFluxGhost::KinetostatFluxGhost(Kinetostat * kinetostat,
 
1922
                                           const string & regulatorPrefix) :
 
1923
    KinetostatFlux(kinetostat,regulatorPrefix)
 
1924
  {
 
1925
    // flag for performing boundary flux calculation
 
1926
    fieldMask_(VELOCITY,FLUX) = false;
 
1927
  }
 
1928
 
 
1929
  //--------------------------------------------------------
 
1930
  //  constructor_transfers
 
1931
  //    instantiates or obtains all dependency managed data
 
1932
  //--------------------------------------------------------
 
1933
  void KinetostatFluxGhost::construct_transfers()
 
1934
  {
 
1935
    KinetostatFlux::construct_transfers();
 
1936
    if (!nodalGhostForce_) {
 
1937
      throw ATC_Error("StressFluxGhost::StressFluxGhost - ghost atoms must be specified");
 
1938
    }
 
1939
  }
 
1940
 
 
1941
  //--------------------------------------------------------
 
1942
  //  compute_boundary_flux:
 
1943
  //    computes the boundary flux to be consistent with
 
1944
  //    the controller
 
1945
  //--------------------------------------------------------
 
1946
  void KinetostatFluxGhost::compute_boundary_flux(FIELDS & fields)
 
1947
  {
 
1948
    // This is only used in computation of atomic sources
 
1949
    boundaryFlux_[VELOCITY] = 0.;
 
1950
  }
 
1951
 
 
1952
  //--------------------------------------------------------
 
1953
  //  add_to_momentum:
 
1954
  //            determines what if any contributions to the
 
1955
  //            finite element equations are needed for
 
1956
  //            consistency with the kinetostat
 
1957
  //--------------------------------------------------------
 
1958
  void KinetostatFluxGhost::add_to_momentum(const DENS_MAT & nodalLambdaForce,
 
1959
                                       DENS_MAT & deltaMomentum,
 
1960
                                       double dt)
 
1961
  {
 
1962
    deltaMomentum.resize(nNodes_,nsd_);
 
1963
    const DENS_MAT & boundaryFlux(nodalGhostForce_->quantity());
 
1964
    for (int i = 0; i < nNodes_; i++) {
 
1965
      for (int j = 0; j < nsd_; j++) {
 
1966
        deltaMomentum(i,j) = nodalLambdaForce(i,j) + dt*boundaryFlux(i,j);
 
1967
      }
 
1968
    }
 
1969
  }
 
1970
 
 
1971
  //--------------------------------------------------------
 
1972
  //  set_kinetostat_rhs
 
1973
  //            sets up the RHS of the kinetostat equations
 
1974
  //            for the coupling parameter lambda
 
1975
  //--------------------------------------------------------
 
1976
  void KinetostatFluxGhost::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
 
1977
  {
 
1978
    // (a) for flux based : 
 
1979
    // form rhs :  \int N_I r dV - \sum_g N_Ig^* f_g
 
1980
    // sources are set in ATC transfer
 
1981
    rhs.reset(nNodes_,nsd_);
 
1982
    const DENS_MAT & momentumSource(momentumSource_.quantity());
 
1983
    const set<int> & applicationNodes(applicationNodes_->quantity());
 
1984
    set<int>::const_iterator iNode;
 
1985
    for (iNode = applicationNodes.begin(); iNode != applicationNodes.end(); iNode++) {
 
1986
      for (int j = 0; j < nsd_; j++) {
 
1987
        rhs(*iNode,j) = momentumSource(*iNode,j);
 
1988
      }
 
1989
    }
 
1990
  }
 
1991
 
 
1992
  //--------------------------------------------------------
 
1993
  //--------------------------------------------------------
 
1994
  //  Class KinetostatFixed
 
1995
  //--------------------------------------------------------
 
1996
  //--------------------------------------------------------
 
1997
 
 
1998
  //--------------------------------------------------------
 
1999
  //  Constructor
 
2000
  //         Grab references to ATC and kinetostat data
 
2001
  //--------------------------------------------------------
 
2002
  KinetostatFixed::KinetostatFixed(Kinetostat * kinetostat,
 
2003
                                   const string & regulatorPrefix) :
 
2004
    KinetostatGlcFs(kinetostat,regulatorPrefix),
 
2005
    mdMassMatrix_(atc_->set_mass_mat_md(VELOCITY)),
 
2006
    isFirstTimestep_(true)
 
2007
  {
 
2008
    // do nothing
 
2009
  }
 
2010
 
 
2011
  //--------------------------------------------------------
 
2012
  //  constructor_transfers
 
2013
  //    instantiates or obtains all dependency managed data
 
2014
  //--------------------------------------------------------
 
2015
  void KinetostatFixed::construct_transfers()
 
2016
  {
 
2017
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
2018
 
 
2019
    // set up node mappings
 
2020
    create_node_maps();
 
2021
 
 
2022
    // determine if map is needed and set up if so
 
2023
    if (this->use_local_shape_functions()) {
 
2024
      lambdaAtomMap_ = new AtomToElementset(atc_,elementMask_);
 
2025
      interscaleManager.add_per_atom_int_quantity(lambdaAtomMap_,
 
2026
                                                  regulatorPrefix_+"LambdaAtomMap");
 
2027
      shapeFunctionMatrix_ = new LocalLambdaCouplingMatrix(atc_,
 
2028
                                                           lambdaAtomMap_,
 
2029
                                                           nodeToOverlapMap_);
 
2030
    }
 
2031
    else {
 
2032
      shapeFunctionMatrix_ = new LambdaCouplingMatrix(atc_,nodeToOverlapMap_);
 
2033
    }
 
2034
    interscaleManager.add_per_atom_sparse_matrix(shapeFunctionMatrix_,
 
2035
                                                 regulatorPrefix_+"LambdaCouplingMatrixMomentum");
 
2036
    linearSolverType_ = AtomicRegulator::CG_SOLVE;
 
2037
 
 
2038
    // base class transfers
 
2039
    KinetostatGlcFs::construct_transfers();
 
2040
  }
 
2041
 
 
2042
  //--------------------------------------------------------
 
2043
  //  initialize
 
2044
  //    initializes all method data
 
2045
  //--------------------------------------------------------
 
2046
  void KinetostatFixed::initialize()
 
2047
  {
 
2048
    KinetostatGlcFs::initialize();
 
2049
 
 
2050
    // reset data to zero
 
2051
    deltaFeMomentum_.reset(nNodes_,nsd_);
 
2052
    deltaNodalAtomicMomentum_.reset(nNodes_,nsd_);
 
2053
  }
 
2054
 
 
2055
  //--------------------------------------------------------
 
2056
  //  halve_force:
 
2057
  //    flag to halve the lambda force for improved
 
2058
  //    accuracy
 
2059
  //--------------------------------------------------------
 
2060
  bool KinetostatFixed::halve_force()
 
2061
  {
 
2062
    if (isFirstTimestep_ || ((atc_->atom_to_element_map_type() == EULERIAN)
 
2063
                             && (atc_->atom_to_element_map_frequency() > 1)
 
2064
                             && (atc_->step() % atc_->atom_to_element_map_frequency() == 1))) {
 
2065
      return true;
 
2066
    }
 
2067
    return false;
 
2068
  }
 
2069
 
 
2070
  //--------------------------------------------------------
 
2071
  //  construct_regulated_nodes:
 
2072
  //    constructs the set of nodes being regulated
 
2073
  //--------------------------------------------------------
 
2074
  void KinetostatFixed::construct_regulated_nodes()
 
2075
  {
 
2076
    InterscaleManager & interscaleManager(atc_->interscale_manager());
 
2077
 
 
2078
    if (!atomicRegulator_->use_localized_lambda()) {
 
2079
      regulatedNodes_ = new RegulatedNodes(atc_);
 
2080
    }
 
2081
    else if (kinetostat_->coupling_mode() == Kinetostat::FLUX) {
 
2082
      regulatedNodes_ = new FixedNodes(atc_);
 
2083
    }
 
2084
    else if (kinetostat_->coupling_mode() == Kinetostat::FIXED) {
 
2085
      // include boundary nodes
 
2086
      regulatedNodes_ = new FixedBoundaryNodes(atc_);
 
2087
    }
 
2088
    else {
 
2089
      throw ATC_Error("ThermostatFixed::construct_regulated_nodes - couldn't determine set of regulated nodes");
 
2090
    }
 
2091
     
 
2092
    interscaleManager.add_set_int(regulatedNodes_,
 
2093
                                  regulatorPrefix_+"RegulatedNodes");
 
2094
 
 
2095
    applicationNodes_ = regulatedNodes_;
 
2096
 
 
2097
    // special set of boundary elements for defining regulated atoms 
 
2098
    if (atomicRegulator_->use_localized_lambda()) {
 
2099
      elementMask_ = new ElementMaskNodeSet(atc_,applicationNodes_);
 
2100
      interscaleManager.add_dense_matrix_bool(elementMask_,
 
2101
                                                   regulatorPrefix_+"BoundaryElementMask");
 
2102
    }
 
2103
  }
 
2104
 
 
2105
  //--------------------------------------------------------
 
2106
  //  initialize_delta_nodal_atomic_momentum:
 
2107
  //    initializes storage for the variable tracking
 
2108
  //    the change in the nodal atomic momentum
 
2109
  //    that has occured over the past timestep
 
2110
  //--------------------------------------------------------
 
2111
  void KinetostatFixed::initialize_delta_nodal_atomic_momentum(double dt)
 
2112
  {
 
2113
    // initialize delta energy
 
2114
    const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity());
 
2115
    initialNodalAtomicMomentum_ = myNodalAtomicMomentum;
 
2116
    initialNodalAtomicMomentum_ *= -1.; // initially stored as negative for efficiency
 
2117
    timeFilter_->apply_pre_step1(nodalAtomicMomentumFiltered_.set_quantity(),
 
2118
                                 myNodalAtomicMomentum,dt);
 
2119
  }
 
2120
 
 
2121
  //--------------------------------------------------------
 
2122
  //  compute_delta_nodal_atomic_momentum:
 
2123
  //    computes the change in the nodal atomic momentum
 
2124
  //    that has occured over the past timestep
 
2125
  //--------------------------------------------------------
 
2126
  void KinetostatFixed::compute_delta_nodal_atomic_momentum(double dt)
 
2127
  {
 
2128
    // set delta energy based on predicted atomic velocities
 
2129
    const DENS_MAT & myNodalAtomicMomentum(nodalAtomicMomentum_->quantity());
 
2130
    timeFilter_->apply_post_step1(nodalAtomicMomentumFiltered_.set_quantity(),
 
2131
                                  myNodalAtomicMomentum,dt);
 
2132
    deltaNodalAtomicMomentum_ = initialNodalAtomicMomentum_;
 
2133
    deltaNodalAtomicMomentum_ += myNodalAtomicMomentum;
 
2134
  }
 
2135
 
 
2136
  //--------------------------------------------------------
 
2137
  //  apply_pre_predictor:
 
2138
  //    apply the kinetostat to the atoms in the
 
2139
  //    pre-predictor integration phase
 
2140
  //--------------------------------------------------------
 
2141
  void KinetostatFixed::apply_pre_predictor(double dt)
 
2142
  {
 
2143
    // initialize values to be track change in finite element energy over the timestep
 
2144
    initialize_delta_nodal_atomic_momentum(dt);
 
2145
    initialFeMomentum_ = -1.*((mdMassMatrix_.quantity())*(velocity_.quantity())); // initially stored as negative for efficiency
 
2146
 
 
2147
    KinetostatGlcFs::apply_pre_predictor(dt);
 
2148
  }
 
2149
 
 
2150
  //--------------------------------------------------------
 
2151
  //  apply_post_corrector:
 
2152
  //    apply the kinetostat to the atoms in the
 
2153
  //    post-corrector integration phase
 
2154
  //--------------------------------------------------------
 
2155
  void KinetostatFixed::apply_post_corrector(double dt)
 
2156
  {  
 
2157
    KinetostatGlcFs::apply_post_corrector(dt);
 
2158
 
 
2159
    // update filtered momentum with lambda force
 
2160
    DENS_MAT & myNodalAtomicLambdaForce(nodalAtomicLambdaForce_->set_quantity());
 
2161
    timeFilter_->apply_post_step2(nodalAtomicMomentumFiltered_.set_quantity(),
 
2162
                                  myNodalAtomicLambdaForce,dt);
 
2163
 
 
2164
    if (halve_force()) {
 
2165
      // Halve lambda force due to fixed temperature constraints
 
2166
      // 1) makes up for poor initial condition
 
2167
      // 2) accounts for possibly large value of lambda when atomic shape function values change
 
2168
      //    from eulerian mapping after more than 1 timestep
 
2169
      //    avoids unstable oscillations arising from 
 
2170
      //    thermostat having to correct for error introduced in lambda changing the 
 
2171
      //    shape function matrices
 
2172
      *lambda_ *= 0.5;
 
2173
    }
 
2174
    
 
2175
    isFirstTimestep_ = false;
 
2176
  }
 
2177
 
 
2178
  //--------------------------------------------------------
 
2179
  //  compute_kinetostat
 
2180
  //            manages the solution and application of the
 
2181
  //            kinetostat equations and variables
 
2182
  //--------------------------------------------------------
 
2183
  void KinetostatFixed::compute_lambda(double dt)
 
2184
  {
 
2185
    // compute predicted changes in nodal atomic momentum
 
2186
    compute_delta_nodal_atomic_momentum(dt);
 
2187
 
 
2188
    // change in finite element momentum
 
2189
    deltaFeMomentum_ = initialFeMomentum_;
 
2190
    deltaFeMomentum_ += (mdMassMatrix_.quantity())*(velocity_.quantity());
 
2191
 
 
2192
    // set up rhs for lambda equation
 
2193
    KinetostatGlcFs::compute_lambda(dt);
 
2194
  }
 
2195
 
 
2196
  //--------------------------------------------------------
 
2197
  //  set_kinetostat_rhs
 
2198
  //            sets up the RHS of the kinetostat equations
 
2199
  //            for the coupling parameter lambda
 
2200
  //--------------------------------------------------------
 
2201
  void KinetostatFixed::set_kinetostat_rhs(DENS_MAT & rhs, double dt)
 
2202
  {
 
2203
    // for essential bcs (fixed nodes) :
 
2204
    // form rhs : (delUpsV - delUps)/dt
 
2205
    const set<int> & regulatedNodes(regulatedNodes_->quantity());
 
2206
    double factor = (1./dt);
 
2207
    for (int i = 0; i < nNodes_; i++) {
 
2208
      if (regulatedNodes.find(i) != regulatedNodes.end()) {
 
2209
        for (int j = 0; j < nsd_; j++) {
 
2210
          rhs(i,j) = factor*(deltaNodalAtomicMomentum_(i,j) - deltaFeMomentum_(i,j));
 
2211
        }
 
2212
      }
 
2213
      else {
 
2214
        for (int j = 0; j < nsd_; j++) {
 
2215
          rhs(i,j) = 0.;
 
2216
        }
 
2217
      }
 
2218
    }
 
2219
  }
 
2220
 
 
2221
  //--------------------------------------------------------
 
2222
  //  add_to_momentum:
 
2223
  //            determines what if any contributions to the
 
2224
  //            finite element equations are needed for
 
2225
  //            consistency with the kinetostat
 
2226
  //--------------------------------------------------------
 
2227
  void KinetostatFixed::add_to_momentum(const DENS_MAT & nodalLambdaForce,
 
2228
                                       DENS_MAT & deltaMomentum,
 
2229
                                       double dt)
 
2230
  {
 
2231
    deltaMomentum.resize(nNodes_,nsd_);
 
2232
    const set<int> & regulatedNodes(regulatedNodes_->quantity());
 
2233
    for (int i = 0; i < nNodes_; i++) {
 
2234
      if (regulatedNodes.find(i) != regulatedNodes.end()) {
 
2235
        for (int j = 0; j < nsd_; j++) {
 
2236
          deltaMomentum(i,j) = 0.;
 
2237
        }
 
2238
      }
 
2239
      else {
 
2240
        for (int j = 0; j < nsd_; j++) {
 
2241
          deltaMomentum(i,j) = nodalLambdaForce(i,j);
 
2242
        }
 
2243
      }
 
2244
    }
 
2245
  }
1234
2246
};