1 #include "PoissonSolver.h"
2 #include "ATC_Coupling.h"
3 #include "FE_Engine.h"
4 #include "PhysicsModel.h"
5 #include "PrescribedDataManager.h"
6 #include "LinearSolver.h"
7 #include <utility>
8 #include <iostream>
9
10 using std::pair;
11
12
13
14 namespace ATC {
15
16 // ====================================================================
17 // PoissonSolver
18 // ====================================================================
PoissonSolver(const FieldName fieldName,const PhysicsModel * physicsModel,const FE_Engine * feEngine,const PrescribedDataManager * prescribedDataMgr,ATC_Coupling * atc,const Array2D<bool> & rhsMask,const int solverType,bool parallel)19 PoissonSolver::PoissonSolver(
20 const FieldName fieldName,
21 const PhysicsModel * physicsModel,
22 const FE_Engine * feEngine,
23 const PrescribedDataManager * prescribedDataMgr,
24 /*const*/ ATC_Coupling * atc,
25 const Array2D<bool> & rhsMask,
26 const int solverType,
27 bool parallel
28 )
29 : atc_(atc),
30 feEngine_(feEngine),
31 prescribedDataMgr_(prescribedDataMgr),
32 physicsModel_(physicsModel),
33 fieldName_(fieldName),
34 rhsMask_(rhsMask),
35 linear_(false),
36 solver_(NULL),
37 solverNL_(NULL),
38 tangent_(NULL),
39 solverType_(solverType),
40 solverTol_(0),
41 solverMaxIter_(0),
42 integrationType_(FULL_DOMAIN),
43 parallel_(parallel)
44 {
45 if (physicsModel_->has_linear_rhs(fieldName)) {
46 linear_ = true;
47 rhsMask_(fieldName,FLUX) = false;
48 }
49
50 else {
51 rhsMask_(fieldName,FLUX) = true;
52 rhsMask_(fieldName,SOURCE) = true;
53 }
54 if (prescribedDataMgr_->has_robin_source(fieldName)) {
55
56
57
58 rhsMask_(fieldName,ROBIN_SOURCE) = true;
59 }
60 }
61 // --------------------------------------------------------------------
~PoissonSolver()62 PoissonSolver::~PoissonSolver()
63 {
64 if (tangent_) delete tangent_;
65 if (solverNL_) delete solverNL_;
66 if (solver_) delete solver_;
67 }
68
69 // --------------------------------------------------------------------
70 // Parser
71 // --------------------------------------------------------------------
72
modify(int narg,char ** arg)73 bool PoissonSolver::modify(int narg, char **arg)
74 {
75 bool match = false;
76 /*! \page man_poisson_solver fix_modify AtC poisson_solver
77 \section syntax
78 fix_modify AtC poisson_solver mesh create <nx> <ny> <nz> <region-id>
79 <f|p> <f|p> <f|p>
80 - nx ny nz = number of elements in x, y, z
81 - region-id = id of region that is to be meshed
82 - f p p = perioidicity flags for x, y, z
83 \section examples
84 <TT> fix_modify AtC poisson_solver mesh create 10 1 1 feRegion p p p </TT>
85 \section description
86 Creates a uniform mesh in a rectangular region
87 \section restrictions
88 creates only uniform rectangular grids in a rectangular region
89 \section related
90 \section default
91 none
92 */
93 int argIdx = 0;
94 if (strcmp(arg[argIdx],"poisson_solver")==0) {
95 argIdx++;
96 if (strcmp(arg[argIdx],"mesh")==0) {
97 argIdx++;
98 // create a FE_Engine
99
100 //feEngine_ = new FE_Engine(this); need alternate constructor?
101 // send args to new engine
102 // arg[0] = "mesh";
103 // arg[1] = "create";
104 // feEngine_->modify(narg,arg);
105
106 }
107 }
108 return match;
109 }
110 // --------------------------------------------------------------------
111 // Initialize
112 // --------------------------------------------------------------------
initialize(void)113 void PoissonSolver::initialize(void)
114 {
115 nNodes_ = feEngine_->num_nodes();
116
117 if (atc_->source_atomic_quadrature(fieldName_))
118 integrationType_ = FULL_DOMAIN_ATOMIC_QUADRATURE_SOURCE;
119
120 // compute penalty for Dirichlet boundaries
121 if (prescribedDataMgr_->none_fixed(fieldName_))
122 throw ATC_Error("Poisson solver needs Dirichlet data");
123
124 const BC_SET & bcs = (prescribedDataMgr_->bcs(fieldName_))[0];
125
126 if (linear_) { // constant rhs
127 if (! solver_ ) {
128 pair<FieldName,FieldName> row_col(fieldName_,fieldName_);
129 Array2D <bool> rhsMask(NUM_FIELDS,NUM_FLUX);
130 rhsMask = false; rhsMask(fieldName_,FLUX) = true;
131 if (prescribedDataMgr_->has_robin_source(fieldName_)) {
132 rhsMask(fieldName_,ROBIN_SOURCE) = true;
133 }
134 // compute stiffness for Poisson solve
135 atc_->compute_rhs_tangent(row_col, rhsMask, atc_->fields(),
136 stiffness_, FULL_DOMAIN, physicsModel_);
137 // create solver
138 solver_ = new LinearSolver(stiffness_,bcs,solverType_,-1,parallel_);
139 }
140 else {
141
142 // re-initialize
143 solver_->initialize(&bcs);
144 }
145 if (solverTol_) solver_->set_tolerance(solverTol_);
146 if (solverMaxIter_) solver_->set_max_iterations(solverMaxIter_);
147
148 }
149 else {
150 // print_mask(rhsMask_);
151 if ( solverNL_ ) delete solverNL_;
152 tangent_ = new PhysicsModelTangentOperator(atc_,physicsModel_, rhsMask_, integrationType_, fieldName_);
153
154 solverNL_ = new NonLinearSolver(tangent_,&bcs,0,parallel_);
155
156 if (solverTol_) solverNL_->set_residual_tolerance(solverTol_);
157 if (solverMaxIter_) solverNL_->set_max_iterations(solverMaxIter_);
158 }
159 }
160
161 // --------------------------------------------------------------------
162 // Solve
163 // --------------------------------------------------------------------
solve(FIELDS & fields,FIELDS & rhs)164 bool PoissonSolver::solve(FIELDS & fields, FIELDS & rhs)
165 {
166 atc_->compute_rhs_vector(rhsMask_, fields, rhs,
167 integrationType_, physicsModel_);
168 CLON_VEC f = column(fields[fieldName_].set_quantity(),0);
169 CLON_VEC r = column(rhs[fieldName_].quantity(),0);
170 bool converged = false;
171 if (linear_) {converged = solver_->solve(f,r);}
172 else {converged = solverNL_->solve(f);}
173
174 if (atc_->source_atomic_quadrature(fieldName_)
175 && LammpsInterface::instance()->atom_charge() ) set_charges(fields);
176 return converged;
177 }
solve(DENS_MAT & field,const DENS_MAT & rhs)178 bool PoissonSolver::solve(DENS_MAT & field, const DENS_MAT & rhs)
179 {
180
181 CLON_VEC f = column(field,0);
182 CLON_VEC r = column(rhs,0);
183 bool converged = false;
184 if (linear_) {converged = solver_->solve(f,r);}
185 else {converged = solverNL_->solve(f);}
186
187 if (atc_->source_atomic_quadrature(fieldName_)
188 && LammpsInterface::instance()->atom_charge() ) set_charges(atc_->fields());
189 return converged;
190 }
191
192 // --------------------------------------------------------------------
193 // set charges on atoms
194 // --------------------------------------------------------------------
set_charges(FIELDS & fields)195 void PoissonSolver::set_charges(FIELDS & fields)
196 {
197 FIELD_MATS sources;
198
199 atc_->compute_sources_at_atoms(rhsMask_, fields, physicsModel_,sources);
200 FIELD_MATS::const_iterator nField = sources.find(fieldName_);
201 if (nField != sources.end()) {
202 const DENS_MAT & electronCharges = nField->second;
203 double * q = LammpsInterface::instance()->atom_charge();
204 int nLocal = atc_->nlocal();
205 if (nLocal > 0) {
206 const Array<int> & i2a = atc_->internal_to_atom_map();
207 for (int i=0; i < nLocal; i++) {
208
209 int atomIdx = i2a(i);
210 q[atomIdx] = -electronCharges(i,0);
211 }
212 }
213 }
214 }
215
216 } // namespace ATC
217