Mechanics of Defect Evolution Library
Loading...
Searching...
No Matches
GalerkinClimbSolver.cpp
Go to the documentation of this file.
1/* This file is part of MoDELib, the Mechanics Of Defects Evolution Library.
2 *
3 *
4 * MoDELib is distributed without any warranty under the
5 * GNU General Public License (GPL) v2 <http://www.gnu.org/licenses/>.
6 */
7
8#ifndef model_GalerkinClimbSolver_cpp_
9#define model_GalerkinClimbSolver_cpp_
10
11#include <deque>
12
13
14#include <ClusterDynamics/ClusterDynamicsParameters.h>
15#include <GalerkinClimbSolver.h>
16#include <ConstrainedVelocitySolve.h>
17#include <TerminalColors.h>
18#include <EqualIteratorRange.h>
19
20#include <DislocationNetwork.h>
21
22namespace model
23{
24
25 template <int dim>
26 GalerkinClimbSolver<dim>::GalerkinClimbSolver(const DislocationNetwork<dim>& DN_in,const ClusterDynamics<dim>* const CD_in) :
27 /* init */ DislocationClimbSolverBase<dim>(DN_in,CD_in)
28 {
29 std::cout<<greenBoldColor<<"Creating GalerkinClimbSolver"<<defaultColor<<std::endl;
30 }
31
32 template <int dim>
33 typename GalerkinClimbSolver<dim>::ForceVectorMatrixType GalerkinClimbSolver<dim>::clusterForceKernel(const int& k,const NetworkLinkType& fieldSegment) const
34 {
35 const Eigen::Matrix<double,1,mSize> deltaC(fieldSegment.quadraturePoint(k).cDD-fieldSegment.quadraturePoint(k).cCD);
36 const double u(fieldSegment.quadraturePoint(k).abscissa);
37 return (ForceVectorMatrixType()<<(1.0-u)*deltaC,u*deltaC).finished();
38 }
39
40 template <int dim>
41 typename GalerkinClimbSolver<dim>::ForceVectorMatrixType GalerkinClimbSolver<dim>::clusterForceVector(const NetworkLinkType& fieldSegment) const
42 {
43 ForceVectorMatrixType F(ForceVectorMatrixType::Zero());
44 const auto bxt(fieldSegment.burgers().cross(fieldSegment.chord()));
45 NetworkLinkType::QuadratureDynamicType::integrate(fieldSegment.quadraturePoints().size(),this,F,&GalerkinClimbSolver<dim>::clusterForceKernel,fieldSegment);
46 F.row(0)*=bxt.dot(fieldSegment.source->climbDirection());
47 F.row(1)*=bxt.dot(fieldSegment.sink->climbDirection());
48 return F;
49 }
50
51 template <int dim>
52 typename GalerkinClimbSolver<dim>::StiffnessMatrixType GalerkinClimbSolver<dim>::clusterStiffnessKernel(const int& k,const NetworkLinkType& fieldSegment,const NetworkLinkType& sourceSegment) const
53 {
54 StiffnessMatrixType temp(StiffnessMatrixType::Zero());
55 if(sourceSegment.grains().size() == 1)
56 {
57 const auto bxt(fieldSegment.burgers().cross(fieldSegment.chord()));
58 const double u(fieldSegment.quadraturePoint(k).abscissa);
59 const Eigen::Matrix<double,2,1> m(( Eigen::Matrix<double,2,1>()<<(1.0-u)*bxt.dot(fieldSegment.source->climbDirection()),
60 /* */ u*bxt.dot(fieldSegment.sink->climbDirection())).finished());
61 const auto sourceConcentratoinMatrices(sourceSegment.concentrationMatrices(fieldSegment.quadraturePoint(k).r,this->CD->cdp));
62 for(int kc=0; kc < mSize; kc++)
63 {
64 temp.template block<2,2>(kc*2,0)+=m*sourceConcentratoinMatrices.row(kc);
65 }
66 }
67 return temp;
68 }
69
70 template <int dim>
71 typename GalerkinClimbSolver<dim>::StiffnessMatrixType GalerkinClimbSolver<dim>::clusterStiffnessMatrix(const NetworkLinkType& fieldSegment,const NetworkLinkType& sourceSegment) const
72 {
73 StiffnessMatrixType K(StiffnessMatrixType::Zero());
74 NetworkLinkType::QuadratureDynamicType::integrate(fieldSegment.quadraturePoints().size(),this,K,&GalerkinClimbSolver<dim>::clusterStiffnessKernel,fieldSegment,sourceSegment);
75 return K;
76 }
77
78 template <int dim>
79 Eigen::VectorXd GalerkinClimbSolver<dim>::getNodeVelocitiesPipe() const
80 {
81 Eigen::VectorXd nodeVelocities(Eigen::VectorXd::Zero(this->DN.networkNodes().size()*dim));
82 return nodeVelocities;
83 }
84
85 template <int dim>
86 void GalerkinClimbSolver<dim>::computeClimbScalarVelocities(const bool& updateStiffness)
87 {
88 std::cout<<", climbSolver "<<std::flush;
89 computeClimbScalarVelocitiesBulk(updateStiffness);
90 }
91
92 // Retained temporarily for an easy source-level comparison with the
93 // optimized assembly below. The complete original file is also backed up
94 // as GalerkingClimbSolver_BK.cpp.
95#if 0
96 template <int dim>
97 void GalerkinClimbSolver<dim>::computeClimbScalarVelocitiesBulk(const bool& updateStiffness)
98 {
99// std::vector<std::vector<Eigen::Triplet<double>>> lhsT(mSize);
100 std::vector<Eigen::VectorXd> Fc(mSize,Eigen::VectorXd::Zero(this->DN.networkNodes().size()));
101 std::vector<Eigen::VectorXd> KKc(mSize,Eigen::VectorXd::Zero(this->DN.networkNodes().size()));
102// bool useLumpedSolver(true);
103
104#ifdef _OPENMP
105 const size_t nThreads(omp_get_max_threads());
106// std::vector<std::vector<std::vector<Eigen::Triplet<double>>>> lhsTV(nThreads,lhsT);
107 std::vector<std::vector<Eigen::VectorXd>> KKcT(nThreads,KKc);
108 std::vector<std::vector<Eigen::VectorXd>> FcT(nThreads,Fc);
109 const EqualIteratorRange<typename DislocationNetwork<dim>::NetworkLinkContainerType::const_iterator> eir(this->DN.networkLinks().begin(),this->DN.networkLinks().end(),nThreads);
110#pragma omp parallel for
111 for(size_t thread=0;thread<eir.size();++thread)
112 {
113 auto& Fc_ref(FcT[thread]);
114 auto& KKc_ref(KKcT[thread]);
115 for(auto fieldLinkIter=eir[thread].first;fieldLinkIter!=eir[thread].second;++fieldLinkIter)
116// for(const auto& fieldLink : this->DN.networkLinks())
117 {// sum line-integral part of displacement field per segment
118 const auto& fieldLink(*fieldLinkIter);
119#else
120 auto& Fc_ref(Fc);
121 auto& KKc_ref(KKc);
122 for(const auto& fieldLink : this->DN.networkLinks())
123 {
124#endif
126// for(size_t thread=0;thread< eir.size();++thread)
127// {
128// for(auto fieldLinkIter=eir[thread].first;fieldLinkIter!=eir[thread].second;++fieldLinkIter)
130// {// sum line-integral part of displacement field per segment
131// const auto& fieldLink(*fieldLinkIter);
132 if( !fieldLink.second.lock()->hasZeroBurgers()
133 && fieldLink.second.lock()->isSessile()
134 && !fieldLink.second.lock()->isBoundarySegment()
135 && !fieldLink.second.lock()->isGrainBoundarySegment()
136 && fieldLink.second.lock()->chordLength()>FLT_EPSILON
137 )
138 {
139 const size_t i0(fieldLink.second.lock()->source->gID());
140 const size_t i1(fieldLink.second.lock()-> sink->gID());
141
142 const ForceVectorMatrixType fc(clusterForceVector(*fieldLink.second.lock()));
143 for(int kc=0; kc<mSize; ++kc)
144 {
145//#ifdef _OPENMP
146//#pragma omp critical
147//#endif
148// {
149 Fc_ref[kc](i0)+=fc(0,kc);
150 Fc_ref[kc](i1)+=fc(1,kc);
151// }
152 }
153
154 for(const auto& sourceLink : this->DN.networkLinks())
155 {// sum line-integral part of displacement field per segment
156 if( !sourceLink.second.lock()->hasZeroBurgers()
157 && !sourceLink.second.lock()->isBoundarySegment()
158 && !sourceLink.second.lock()->isGrainBoundarySegment()
159 && sourceLink.second.lock()->chordLength()>FLT_EPSILON
160 )
161 {
162 const size_t j0(sourceLink.second.lock()->source->gID());
163 const size_t j1(sourceLink.second.lock()-> sink->gID());
164 const StiffnessMatrixType kcc(clusterStiffnessMatrix(*fieldLink.second.lock(),*sourceLink.second.lock()));
165
166 for(int kc=0; kc<mSize; ++kc)
167 {
168 const Eigen::Matrix<double,2,2> kccs(kcc.template block<2,2>(2*kc,0));
169
170// if(useLumpedSolver)
171// {
172//#ifdef _OPENMP
173//#pragma omp critical
174//#endif
175// {
176 KKc_ref[kc](i0)+=0.5*kccs(0,0)+0.5*kccs(0,1);
177 KKc_ref[kc](j0)+=0.5*kccs(0,0)+0.5*kccs(1,0);
178 KKc_ref[kc](i1)+=0.5*kccs(1,0)+0.5*kccs(1,1);
179 KKc_ref[kc](j1)+=0.5*kccs(0,1)+0.5*kccs(1,1);
180
181
182 //lhsT[kc].emplace_back(i0,i0,0.5*kccs(0,0));
183 //lhsT[kc].emplace_back(j0,j0,0.5*kccs(0,0));
184
185 //lhsT[kc].emplace_back(i0,i0,0.5*kccs(0,1));
186// lhsT[kc].emplace_back(j1,j1,0.5*kccs(0,1));
187
188// lhsT[kc].emplace_back(i1,i1,0.5*kccs(1,0));
189 //lhsT[kc].emplace_back(j0,j0,0.5*kccs(1,0));
190
191 // lhsT[kc].emplace_back(i1,i1,0.5*kccs(1,1));
192// lhsT[kc].emplace_back(j1,j1,0.5*kccs(1,1));
193// }
194// }
195// else
196// {
197// const bool forceSym(true);// at the moment we need to symmtrize, since numerical intgration on fieldLink is not equivalent to analytical integration on sourceLink
198// if(forceSym)
199// {
200//#ifdef _OPENMP
201//#pragma omp critical
202//#endif
203// {
204// lhsT[kc].emplace_back(i0,j0,0.5*kccs(0,0));
205// lhsT[kc].emplace_back(i0,j1,0.5*kccs(0,1));
206// lhsT[kc].emplace_back(i1,j0,0.5*kccs(1,0));
207// lhsT[kc].emplace_back(i1,j1,0.5*kccs(1,1));
208//
209// lhsT[kc].emplace_back(j0,i0,0.5*kccs(0,0));
210// lhsT[kc].emplace_back(j1,i0,0.5*kccs(0,1));
211// lhsT[kc].emplace_back(j0,i1,0.5*kccs(1,0));
212// lhsT[kc].emplace_back(j1,i1,0.5*kccs(1,1));
213// }
214// }
215// else
216// {
217//#ifdef _OPENMP
218//#pragma omp critical
219//#endif
220// {
221// lhsT[kc].emplace_back(i0,j0,kccs(0,0));
222// lhsT[kc].emplace_back(i0,j1,kccs(0,1));
223// lhsT[kc].emplace_back(i1,j0,kccs(1,0));
224// lhsT[kc].emplace_back(i1,j1,kccs(1,1));
225// }
226// }
227// }
228 }
229 }
230 }
231 }
232 }
233#ifdef _OPENMP
234 }
235 for(size_t thread=0;thread< eir.size();++thread)
236 {// recombine contributions of different threads
237 for(int kc=0; kc<mSize; ++kc)
238 {
239 Fc[kc]+=FcT[thread][kc];
240 KKc[kc]+=KKcT[thread][kc];
241// lhsT[kc].insert(lhsT[kc].end(), lhsTV[thread][kc].begin(), lhsTV[thread][kc].end());
242 }
243 }
244#endif
245
246// Eigen::SparseMatrix<double> Kcc(this->DN.networkNodes().size(),this->DN.networkNodes().size());
247 std::vector<Eigen::Array<double,1,mSize>> nodeV(this->DN.networkNodes().size(),Eigen::Array<double,1,mSize>::Zero());
248
249 for(int kc=0; kc<mSize; ++kc)
250 {
251// Kcc.setFromTriplets(lhsT[kc].begin(),lhsT[kc].end());
252
253// if(size_t(Kcc.rows())!=this->DN.networkNodes().size() || size_t(Kcc.cols())!=this->DN.networkNodes().size())
254// {
255// throw std::runtime_error("the Stiffness Matrix size is not equal to the node size.");
256// }
257
258// if(useLumpedSolver)
259// {
260// Eigen::VectorXd Kccd(Kcc.diagonal());
261 for (size_t n=0; n<this->DN.networkNodes().size(); n++)
262 {
263 if(std::fabs(KKc[kc](n))>FLT_EPSILON)
264 {
265 nodeV[n](kc)=Fc[kc](n)/KKc[kc](n);
266 }
267 }
268// }
269// else
270// {
271// throw std::runtime_error("GalerkinClimbSolver: only the lumped solver is implemented.");
272// }
273 }
274
275 this->scalarVelocities().resize(this->DN.networkNodes().size(),Eigen::Array<double,1,mSize>::Zero());
276 for(size_t k=0; k<this->DN.networkNodes().size();++k)
277 {
278 this->scalarVelocities()[k]=nodeV[k];
279 }
280 }
281
282#endif
283
284 template <int dim>
285 void GalerkinClimbSolver<dim>::computeClimbScalarVelocitiesBulk(const bool& updateStiffness)
286 {
287 const size_t nodeCount(this->DN.networkNodes().size());
288 Fc.assign(mSize,Eigen::VectorXd::Zero(nodeCount));
289 if(updateStiffness)
290 {
291 KKc.assign(mSize,Eigen::VectorXd::Zero(nodeCount));
292 }
293
294 using NetworkLinkPtrType=std::shared_ptr<NetworkLinkType>;
295 struct SourceData
296 {
297 NetworkLinkPtrType link;
298 size_t j0;
299 size_t j1;
300 };
301 struct FieldQuadratureData
302 {
303 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
304 VectorDim r;
305 double weightedM0;
306 double weightedM1;
307 };
308 using FieldQuadratureContainer=std::vector<FieldQuadratureData,Eigen::aligned_allocator<FieldQuadratureData>>;
309 struct FieldData
310 {
311 NetworkLinkPtrType link;
312 size_t i0;
313 size_t i1;
314 ForceVectorMatrixType force;
315 FieldQuadratureContainer quadrature;
316 };
317
318 // Lock and filter each network link once, outside the O(N^2) loop.
319 std::vector<SourceData> sources;
320 std::vector<FieldData> fields;
321 sources.reserve(this->DN.networkLinks().size());
322 fields.reserve(this->DN.networkLinks().size());
323 for(const auto& networkLink : this->DN.networkLinks())
324 {
325 const auto link(networkLink.second.lock());
326 if(!link)
327 {
328 continue;
329 }
330
331 const bool isActive(!link->hasZeroBurgers()
332 && !link->isBoundarySegment()
333 && !link->isGrainBoundarySegment()
334 && link->chordLength()>FLT_EPSILON);
335 // clusterStiffnessKernel returned zero for all other sources.
336 if(isActive && link->grains().size()==1)
337 {
338 sources.push_back({link,link->source->gID(),link->sink->gID()});
339 }
340
341 if(isActive && link->isSessile())
342 {
343 FieldData field{link,link->source->gID(),link->sink->gID(),clusterForceVector(*link),{}};
344 field.quadrature.reserve(link->quadraturePoints().size());
345
346 // These quantities depend only on the field segment and its
347 // quadrature point, so compute them once rather than per source.
348 const auto bxt(link->burgers().cross(link->chord()));
349 const double sourceProjection(bxt.dot(link->source->climbDirection()));
350 const double sinkProjection(bxt.dot(link->sink->climbDirection()));
351 for(const auto& qp : link->quadraturePoints())
352 {
353 field.quadrature.push_back({qp.r,
354 qp.weight*(1.0-qp.abscissa)*sourceProjection,
355 qp.weight*qp.abscissa*sinkProjection});
356 }
357 fields.emplace_back(std::move(field));
358 }
359 }
360
361#ifdef _OPENMP
362 const size_t nThreads(omp_get_max_threads());
363 std::vector<std::vector<Eigen::VectorXd>> KKcT(nThreads,KKc);
364 std::vector<std::vector<Eigen::VectorXd>> FcT(nThreads,Fc);
365#pragma omp parallel for schedule(guided)
366 for(size_t fieldIndex=0;fieldIndex<fields.size();++fieldIndex)
367 {
368 const size_t thread(omp_get_thread_num());
369 auto& Fc_ref(FcT[thread]);
370 auto& KKc_ref(KKcT[thread]);
371 const auto& field(fields[fieldIndex]);
372#else
373 auto& Fc_ref(Fc);
374 auto& KKc_ref(KKc);
375 for(const auto& field : fields)
376 {
377#endif
378 for(int kc=0;kc<mSize;++kc)
379 {
380 Fc_ref[kc](field.i0)+=field.force(0,kc);
381 Fc_ref[kc](field.i1)+=field.force(1,kc);
382 }
383
384 if(updateStiffness)
385 {
386 for(const auto& source : sources)
387 {
388 Eigen::Array<double,mSize,1> i0Term(Eigen::Array<double,mSize,1>::Zero());
389 Eigen::Array<double,mSize,1> j0Term(Eigen::Array<double,mSize,1>::Zero());
390 Eigen::Array<double,mSize,1> i1Term(Eigen::Array<double,mSize,1>::Zero());
391 Eigen::Array<double,mSize,1> j1Term(Eigen::Array<double,mSize,1>::Zero());
392
393 // Assemble exactly the four lumped contributions consumed by
394 // the solver; do not materialize the 2*mSize-by-2 pair matrix.
395 for(const auto& qp : field.quadrature)
396 {
397 const auto concentration(source.link->concentrationMatrices(qp.r,this->CD->cdp));
398 const auto c0(concentration.col(0).array());
399 const auto c1(concentration.col(1).array());
400 i0Term+=0.5*qp.weightedM0*(c0+c1);
401 j0Term+=0.5*(qp.weightedM0+qp.weightedM1)*c0;
402 i1Term+=0.5*qp.weightedM1*(c0+c1);
403 j1Term+=0.5*(qp.weightedM0+qp.weightedM1)*c1;
404 }
405
406 for(int kc=0;kc<mSize;++kc)
407 {
408 KKc_ref[kc](field.i0)+=i0Term(kc);
409 KKc_ref[kc](source.j0)+=j0Term(kc);
410 KKc_ref[kc](field.i1)+=i1Term(kc);
411 KKc_ref[kc](source.j1)+=j1Term(kc);
412 }
413 }
414 }
415 }
416
417#ifdef _OPENMP
418 for(size_t thread=0;thread<nThreads;++thread)
419 {
420 for(int kc=0;kc<mSize;++kc)
421 {
422 Fc[kc]+=FcT[thread][kc];
423 if(updateStiffness)
424 {
425 KKc[kc]+=KKcT[thread][kc];
426 }
427 }
428 }
429#endif
430
431 const auto Z=periodicClimbVelocityMap<dim>(this->DN);
432 std::vector<Eigen::Array<double,1,mSize>> nodeV(nodeCount,Eigen::Array<double,1,mSize>::Zero());
433 for(int kc=0;kc<mSize;++kc)
434 {
435 const Eigen::VectorXd constrained=solveConstrainedVelocityDiagonal(KKc[kc],Fc[kc],Z);
436 for(size_t n=0;n<nodeCount;++n)
437 {
438 nodeV[n](kc)=constrained(n);
439 }
440 }
441
442 this->scalarVelocities().resize(nodeCount,Eigen::Array<double,1,mSize>::Zero());
443 for(size_t k=0;k<nodeCount;++k)
444 {
445 this->scalarVelocities()[k]=nodeV[k];
446 }
447 }
448
449 template <int dim>
450 Eigen::VectorXd GalerkinClimbSolver<dim>::getNodeVelocitiesBulk() const
451 {
452 Eigen::VectorXd nodeVelocities(Eigen::VectorXd::Zero(this->DN.networkNodes().size()*dim));
453 size_t k=0;
454 for(auto& node: this->DN.networkNodes())
455 {
456 const double vScalarTotal(-1.0*(this->scalarVelocities()[k]*this->CD->cdp.msVector/this->CD->cdp.msVector.abs()).matrix().sum());
457 nodeVelocities.template segment<dim>(k*dim)= vScalarTotal*node.second.lock()->climbDirection();
458 k++;
459 }
460 return interpolateClimbSubdivisionVelocities<dim>(this->DN,nodeVelocities);
461 }
462
463 template <int dim>
464 Eigen::VectorXd GalerkinClimbSolver<dim>::getNodeVelocities() const
465 {
466 return getNodeVelocitiesBulk()+getNodeVelocitiesPipe();
467 }
468
469 template class GalerkinClimbSolver<3>;
470
471}
472#endif
model::GalerkinClimbSolver::GalerkinClimbSolver
GalerkinClimbSolver(const DislocationNetwork< dim > &, const ClusterDynamics< dim > *const)
create_polycrystal_figure.source
str source
Definition create_polycrystal_figure.py:11
model::controls::field
void field(QGridLayout *layout, int row, const QString &title, QWidget *input, const QString &description)
Definition DDqtControls.h:51
model::solveConstrainedVelocityDiagonal
Eigen::VectorXd solveConstrainedVelocityDiagonal(const Eigen::VectorXd &diagonal, const Eigen::VectorXd &force, const VelocityMatrix &Z)
Definition ConstrainedVelocitySolve.h:15
plotLoops.r
float r
Definition plotLoops.py:41
VectorDim
Eigen::Matrix< double, 3, 1 > VectorDim
Definition pyMoDELib.cxx:66
u
REAL u[4]
Definition triangle.cpp:5168
k
int k
Definition triangle.cpp:3300
temp
vertex temp
Definition triangle.cpp:9214
first
vertex first
Definition triangle.cpp:8231
m
void t struct mesh * m
Definition triangle.cpp:3669
Generated by doxygen 1.15.0