← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/atc/ImplicitSolveOperator.cpp

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/atc/ImplicitSolveOperator.cpp
7
7
#include "PhysicsModel.h"
8
8
#include "PrescribedDataManager.h"
9
9
 
 
10
 
 
11
using std::set;
 
12
using std::vector;
 
13
 
10
14
namespace ATC {
11
15
 
12
16
// --------------------------------------------------------------------
15
19
// --------------------------------------------------------------------
16
20
// --------------------------------------------------------------------
17
21
ImplicitSolveOperator::
18
 
ImplicitSolveOperator(ATC_Coupling * atc,
19
 
                      /*const*/ FE_Engine * feEngine,
20
 
                      const PhysicsModel * physicsModel)
21
 
  : atc_(atc),
22
 
    feEngine_(feEngine),
23
 
    physicsModel_(physicsModel)
 
22
ImplicitSolveOperator(double alpha, double dt)
 
23
  : n_(0),
 
24
    dof_(0),
 
25
    dt_(dt),
 
26
    alpha_(alpha),
 
27
    epsilon0_(1.0e-8)
24
28
{
25
29
  // Nothing else to do here
26
30
}
27
31
 
28
32
// --------------------------------------------------------------------
 
33
//  operator *
 
34
// --------------------------------------------------------------------
 
35
DENS_VEC
 
36
ImplicitSolveOperator::operator * (const DENS_VEC & x) const
 
37
{
 
38
  // This method uses a matrix-free approach to approximate the
 
39
  // multiplication by matrix A in the matrix equation Ax=b, where the
 
40
  // matrix equation results from an implicit treatment of the
 
41
  // fast field. In brief, if the ODE for the fast field can be written:
 
42
  //
 
43
  //  dx/dt = R(x)
 
44
  // 
 
45
  // A generalized discretization can be written:
 
46
  //
 
47
  //  1/dt * (x^n+1 - x^n) = alpha * R(x^n+1) + (1-alpha) * R(x^n)
 
48
  //
 
49
  // Taylor expanding the R(x^n+1) term and rearranging gives the
 
50
  // equation to be solved for dx at each timestep:
 
51
  //
 
52
  //  [1 - dt * alpha * dR/dx] * dx = dt * R(x^n)
 
53
  //
 
54
  // The operator defined in this method computes the left-hand side,
 
55
  // given a vector dx.  It uses a finite difference, matrix-free
 
56
  // approximation of dR/dx * dx, giving:
 
57
  //
 
58
  //  [1 - dt * alpha * dR/dx] * dx = dt * R(x^n)
 
59
  //      ~=  dx - dt*alpha/epsilon * ( R(x^n + epsilon*dx) - R(x^n) )
 
60
  //
 
61
  // Compute epsilon
 
62
  double epsilon = (x.norm()>0.0) ? epsilon0_*x0_.norm()/x.norm():epsilon0_;
 
63
  // Compute incremented vector x^n+1 = x^n + epsilon*dx
 
64
  x_ = x0_ + epsilon * x;
 
65
  // Evaluate R(x)
 
66
  this->R(x_,R_);
 
67
  // Compute full left hand side and return it
 
68
  DENS_VEC Ax = x - dt_ * alpha_ / epsilon * (R_ - R0_);
 
69
  return Ax;
 
70
}
 
71
 
 
72
// --------------------------------------------------------------------
 
73
//  rhs of Ax = r
 
74
// --------------------------------------------------------------------
 
75
DENS_VEC
 
76
ImplicitSolveOperator::r() const
 
77
{
 
78
  return dt_ * R0_; // dt * R(T^n)
 
79
}
 
80
 
 
81
// --------------------------------------------------------------------
 
82
//  preconditioner
 
83
// --------------------------------------------------------------------
 
84
DIAG_MAT
 
85
ImplicitSolveOperator::preconditioner() const
 
86
{
 
87
  DENS_VEC diag(n_);
 
88
  diag = 1.0;
 
89
  DIAG_MAT preconditioner(diag);
 
90
  return preconditioner;
 
91
}
 
92
 
 
93
// --------------------------------------------------------------------
29
94
// --------------------------------------------------------------------
30
95
//  FieldImplicitSolveOperator
31
96
// --------------------------------------------------------------------
32
97
// --------------------------------------------------------------------
33
98
FieldImplicitSolveOperator::
34
99
FieldImplicitSolveOperator(ATC_Coupling * atc,
35
 
                           /*const*/ FE_Engine * feEngine,
36
100
                           FIELDS & fields,
37
101
                           const FieldName fieldName,
38
102
                           const Array2D< bool > & rhsMask,
40
104
                           double simTime,
41
105
                           double dt,
42
106
                           double alpha)
43
 
  : ImplicitSolveOperator(atc, feEngine, physicsModel),
 
107
  : ImplicitSolveOperator(alpha, dt),
44
108
    fieldName_(fieldName),
45
 
    fields_(fields), // ref to fields
46
 
    time_(simTime),
47
 
    dt_(dt),
48
 
    alpha_(alpha),
49
 
    epsilon0_(1.0e-8)
 
109
    atc_(atc),
 
110
    physicsModel_(physicsModel),
 
111
    fields0_(fields), // ref to fields
 
112
    fields_ (fields), // copy of fields
 
113
    rhsMask_(rhsMask),
 
114
    time_(simTime)
50
115
{
51
 
  // find field associated with ODE
 
116
  const DENS_MAT & f = fields0_[fieldName_].quantity();
 
117
  dof_ = f.nCols();
 
118
  if (dof_ > 1) throw ATC_Error("Implicit solver operator can only handle scalar fields");
 
119
  // create all to free map
 
120
  int nNodes = f.nRows();
 
121
  set<int> fixedNodes_ = atc_->prescribed_data_manager()->fixed_nodes(fieldName_);
 
122
  n_ = nNodes;
 
123
  vector<bool> tag(nNodes); 
 
124
  set<int>::const_iterator it;  int i = 0;
 
125
  for (i = 0; i < nNodes; ++i) { tag[i] = true; }
 
126
  for (it=fixedNodes_.begin();it!=fixedNodes_.end();++it) {tag[*it]=false;}
 
127
  int m = 0;
 
128
  for (i = 0; i < nNodes; ++i) { if (tag[i]) freeNodes_[i]= m++; }
 
129
//std::cout << " nodes " << n_ << " " << nNodes << "\n";
 
130
  
 
131
  // Save current field
 
132
  x0_.reset(n_);
 
133
  to_free(f,x0_);
 
134
  x_  = x0_; // initialize
 
135
 
 
136
  // righthand side/forcing vector
52
137
  rhsMask_.reset(NUM_FIELDS,NUM_FLUX);
53
138
  rhsMask_ = false;
54
139
  for (int i = 0; i < rhsMask.nCols(); i++) {
55
140
    rhsMask_(fieldName_,i) =  rhsMask(fieldName_,i);
56
141
  }
 
142
//std::cout << print_mask(rhsMask_) << "\n";
57
143
  massMask_.reset(1);
58
144
  massMask_(0) = fieldName_;
59
 
 
60
 
  // Save off current field
61
 
  TnVect_ = column(fields_[fieldName_].quantity(),0); 
62
 
 
63
 
  // Allocate vectors for fields and rhs
64
 
  int nNodes = atc_->num_nodes();
65
 
  // copy fields
66
 
  fieldsNp1_ = fields_;
67
 
  // size rhs
68
 
  int dof = fields_[fieldName_].nCols();
69
 
  RnMap_ [fieldName_].reset(nNodes,dof);
70
 
  RnpMap_[fieldName_].reset(nNodes,dof);
71
 
 
 
145
  rhs_[fieldName_].reset(nNodes,dof_);
72
146
  // Compute the RHS vector R(T^n) 
73
 
  // Set BCs on Rn, multiply by inverse mass and then extract its vector
74
 
  atc_->compute_rhs_vector(rhsMask_, fields_, RnMap_,
75
 
                                   FULL_DOMAIN, physicsModel_);
76
 
  DENS_MAT & Rn = RnMap_[fieldName_].set_quantity();
77
 
  atc_->prescribed_data_manager()
78
 
    ->set_fixed_dfield(time_, fieldName_, Rn);
79
 
  atc_->apply_inverse_mass_matrix(Rn,fieldName_);
80
 
  RnVect_ = column(Rn,0);
81
 
}
82
 
 
83
 
 
84
 
// --------------------------------------------------------------------
85
 
//  operator *(Vector)
86
 
// --------------------------------------------------------------------
87
 
DENS_VEC
88
 
FieldImplicitSolveOperator::operator * (DENS_VEC x) const
89
 
{
90
 
  // This method uses a matrix-free approach to approximate the
91
 
  // multiplication by matrix A in the matrix equation Ax=b, where the
92
 
  // matrix equation results from an implicit treatment of the
93
 
  // fast field solve for the Two Temperature Model.  In
94
 
  // brief, if the ODE for the fast field can be written:
95
 
  //
96
 
  //  dT/dt = R(T)
97
 
  // 
98
 
  // A generalized discretization can be written:
99
 
  //
100
 
  //  1/dt * (T^n+1 - T^n) = alpha * R(T^n+1) + (1-alpha) * R(T^n)
101
 
  //
102
 
  // Taylor expanding the R(T^n+1) term and rearranging gives the
103
 
  // equation to be solved for dT at each timestep:
104
 
  //
105
 
  //  [1 - dt * alpha * dR/dT] * dT = dt * R(T^n)
106
 
  //
107
 
  // The operator defined in this method computes the left-hand side,
108
 
  // given a vector dT.  It uses a finite difference, matrix-free
109
 
  // approximation of dR/dT * dT, giving:
110
 
  //
111
 
  //  [1 - dt * alpha * dR/dT] * dT = dt * R(T^n)
112
 
  //      ~=  dT - dt*alpha/epsilon * ( R(T^n + epsilon*dT) - R(T^n) )
113
 
  //
114
 
  
115
 
  
116
 
  // Compute epsilon
117
 
  double epsilon = (x.norm() > 0.0) ? epsilon0_ * TnVect_.norm()/x.norm() : epsilon0_;
118
 
 
119
 
  // Compute incremented vector = T + epsilon*dT
120
 
  fieldsNp1_[fieldName_] = TnVect_ + epsilon * x;
121
 
 
122
 
  // Evaluate R(b)
123
 
  atc_->compute_rhs_vector(rhsMask_, fieldsNp1_, RnpMap_,
124
 
                                   FULL_DOMAIN, physicsModel_);
125
 
  DENS_MAT & Rnp = RnpMap_[fieldName_].set_quantity();
126
 
  atc_->prescribed_data_manager()
127
 
    ->set_fixed_dfield(time_, fieldName_, Rnp);
128
 
  atc_->apply_inverse_mass_matrix(Rnp,fieldName_);
129
 
  RnpVect_ = column(Rnp,0);
130
 
 
131
 
  // Compute full left hand side and return it
132
 
  DENS_VEC Ax = x - dt_ * alpha_ / epsilon * (RnpVect_ - RnVect_);
133
 
  return Ax;
134
 
}
135
 
 
136
 
// --------------------------------------------------------------------
137
 
//  rhs
138
 
// --------------------------------------------------------------------
139
 
DENS_VEC
140
 
FieldImplicitSolveOperator::rhs()
141
 
{
142
 
  // Return dt * R(T^n)
143
 
  return dt_ * RnVect_;
144
 
}
145
 
 
146
 
// --------------------------------------------------------------------
147
 
//  preconditioner
148
 
// --------------------------------------------------------------------
149
 
DIAG_MAT
150
 
FieldImplicitSolveOperator::preconditioner(FIELDS & fields)
151
 
{
152
 
  // Just create and return identity matrix
153
 
  int nNodes = atc_->num_nodes();
154
 
  DENS_VEC ones(nNodes);
155
 
  ones = 1.0;
156
 
  DIAG_MAT identity(ones);
157
 
  return identity;
158
 
}
 
147
  R0_.reset(n_);
 
148
  R_ .reset(n_);
 
149
  R(x0_, R0_);
 
150
}
 
151
 
 
152
void FieldImplicitSolveOperator::to_all(const VECTOR &x, MATRIX &f) const
 
153
{
 
154
  f.reset(x.nRows(),1);
 
155
  for (int i = 0; i < x.nRows(); ++i) {
 
156
    f(i,0) = x(i);
 
157
  }
 
158
}
 
159
void FieldImplicitSolveOperator::to_free(const MATRIX &r, VECTOR &v) const
 
160
{
 
161
  v.reset(r.nRows());
 
162
  for (int i = 0; i < r.nRows(); ++i) {
 
163
    v(i) = r(i,0);
 
164
  }
 
165
}
 
166
void 
 
167
FieldImplicitSolveOperator::R(const DENS_VEC &x, DENS_VEC &v ) const
 
168
{
 
169
  DENS_MAT & f = fields_[fieldName_].set_quantity();
 
170
  atc_->prescribed_data_manager()->set_fixed_field(time_, fieldName_, f);
 
171
  to_all(x,f);
 
172
  atc_->compute_rhs_vector(rhsMask_,fields_,rhs_,FULL_DOMAIN,physicsModel_);
 
173
  DENS_MAT & r = rhs_[fieldName_].set_quantity();
 
174
  atc_->prescribed_data_manager()->set_fixed_dfield(time_, fieldName_, r);
 
175
  atc_->apply_inverse_mass_matrix(r,fieldName_);
 
176
  to_free(r,v);
 
177
#if 0
 
178
int n = 6;
 
179
//std::cout << "# x "; for (int i = 0; i < n_; ++i)  std::cout << x(i) << " "; std::cout << "\n";
 
180
//std::cout << "# f "; for (int i = 0; i < n; ++i)  std::cout << f(i,0) << " "; std::cout << "\n";
 
181
std::cout << "# r "; for (int i = 0; i < n; ++i)  std::cout << r(i,0) << " "; std::cout << "\n";
 
182
//std::cout << "# v "; for (int i = 0; i < n; ++i)  std::cout << v(i) << " "; std::cout << "\n";
 
183
#endif
 
184
}
 
185
 
 
186
void FieldImplicitSolveOperator::solution(const DENS_MAT & dx, DENS_MAT &f) const
 
187
{
 
188
  DENS_MAT & df = fields_[fieldName_].set_quantity();
 
189
  to_all(column(dx,0),df);
 
190
  atc_->prescribed_data_manager()->set_fixed_dfield(time_, fieldName_, df);
 
191
  f += df;
 
192
}
 
193
void FieldImplicitSolveOperator::rhs(const DENS_MAT & r, DENS_MAT &rhs) const
 
194
{
 
195
  
 
196
  to_all(column(r,0),rhs);
 
197
}
 
198
 
159
199
 
160
200
} // namespace ATC