/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Contents of /trunk/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6799 - (show annotations)
Mon Mar 25 05:53:58 2019 UTC (3 weeks, 4 days ago) by aellery
File size: 65726 byte(s)
I have rewritten the solverbuddy. Briefly:

1. The remaining AMG code has been removed from PASO.
2. If Trilinos is available, escript will now use it by default.
3. eScript will use a direct solver by default, (if one is available,) when solving 2 dimensional problems and an iterative solver, by default, when solving 3 dimensional problems. This can be changed by a user by manually specifying which solver to use.
4. There is a new option available, setHermitian(), that allows a user to specify when a coefficient matrix is Hermitian.
5. Symmetry information is always passed onto the Trilinos solver when this information is relevant.
6. All tests have been updated, when relevant, to reflect these changes.
7. I fixed a couple of undocumented bugs.


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2018 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #include <ripley/RipleyDomain.h>
18 #include <ripley/domainhelpers.h>
19
20 #include <escript/DataFactory.h>
21 #include <escript/FunctionSpaceFactory.h>
22 #include <escript/index.h>
23 #include <escript/SolverOptions.h>
24
25 #ifdef ESYS_HAVE_TRILINOS
26 #include <trilinoswrap/TrilinosMatrixAdapter.h>
27 #endif
28
29 #ifdef ESYS_HAVE_PASO
30 #include <paso/SystemMatrix.h>
31 #include <paso/Transport.h>
32 #endif
33
34 #include <iomanip>
35 #include <iostream>
36
37 namespace bp = boost::python;
38
39 using namespace std;
40 using escript::ValueError;
41 using escript::NotImplementedError;
42
43 #ifdef ESYS_HAVE_TRILINOS
44 using esys_trilinos::TrilinosMatrixAdapter;
45 using esys_trilinos::const_TrilinosGraph_ptr;
46 #endif
47
48 namespace ripley {
49
50 DecompositionPolicy RipleyDomain::m_decompPolicy = DECOMP_ADD_ELEMENTS;
51
52 void RipleyDomain::setDecompositionPolicy(DecompositionPolicy value)
53 {
54 m_decompPolicy = value;
55 }
56
57 DecompositionPolicy RipleyDomain::getDecompositionPolicy()
58 {
59 return m_decompPolicy;
60 }
61
62 void tupleListToMap(DataMap& mapping, const bp::list& list)
63 {
64 for (int i = 0; i < len(list); i++) {
65 if (!bp::extract<bp::tuple>(list[i]).check())
66 throw ValueError("Passed in list contains objects"
67 " other than tuples");
68 bp::tuple t = bp::extract<bp::tuple>(list[i]);
69 if (len(t) != 2 || !bp::extract<string>(t[0]).check() ||
70 !bp::extract<escript::Data>(t[1]).check())
71 throw ValueError("The passed in list must contain tuples"
72 " of the form (string, escript::Data)");
73 mapping[bp::extract<string>(t[0])] = bp::extract<escript::Data>(t[1]);
74 }
75 }
76
77 RipleyDomain::RipleyDomain(dim_t dim, escript::SubWorld_ptr p) :
78 m_numDim(dim),
79 m_status(0)
80 {
81 if (p.get() == NULL)
82 m_mpiInfo = escript::makeInfo(MPI_COMM_WORLD);
83 else
84 m_mpiInfo = p->getMPI();
85
86 assembler_type = DEFAULT_ASSEMBLER;
87 }
88
89 RipleyDomain::~RipleyDomain()
90 {
91 // cleanup of MPI is dealt with by shared_ptr
92 }
93
94 bool RipleyDomain::operator==(const AbstractDomain& other) const
95 {
96 const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
97 if (o) {
98 return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
99 && m_elementTags==o->m_elementTags
100 && m_faceTags==o->m_faceTags);
101 }
102 return false;
103 }
104
105 bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
106 {
107 switch (fsType) {
108 case DegreesOfFreedom:
109 case ReducedDegreesOfFreedom:
110 case Nodes:
111 case ReducedNodes:
112 case Elements:
113 case ReducedElements:
114 case FaceElements:
115 case ReducedFaceElements:
116 case Points:
117 return true;
118 default:
119 break;
120 }
121 return false;
122 }
123
124 string RipleyDomain::functionSpaceTypeAsString(int fsType) const
125 {
126 switch (fsType) {
127 case DegreesOfFreedom:
128 return "Ripley_DegreesOfFreedom [Solution(domain)]";
129 case ReducedDegreesOfFreedom:
130 return "Ripley_ReducedDegreesOfFreedom [ReducedSolution(domain)]";
131 case Nodes:
132 return "Ripley_Nodes [ContinuousFunction(domain)]";
133 case ReducedNodes:
134 return "Ripley_ReducedNodes [ReducedContinuousFunction(domain)]";
135 case Elements:
136 return "Ripley_Elements [Function(domain)]";
137 case ReducedElements:
138 return "Ripley_ReducedElements [ReducedFunction(domain)]";
139 case FaceElements:
140 return "Ripley_FaceElements [FunctionOnBoundary(domain)]";
141 case ReducedFaceElements:
142 return "Ripley_ReducedFaceElements [ReducedFunctionOnBoundary(domain)]";
143 case Points:
144 return "Ripley_Points [DiracDeltaFunctions(domain)]";
145 default:
146 break;
147 }
148 return "Invalid function space type code";
149 }
150
151 pair<int,dim_t> RipleyDomain::getDataShape(int fsType) const
152 {
153 const int ptsPerSample = (m_numDim==2 ? 4 : 8);
154 switch (fsType) {
155 case Nodes:
156 case ReducedNodes:
157 return pair<int,dim_t>(1, getNumNodes());
158 case DegreesOfFreedom:
159 case ReducedDegreesOfFreedom:
160 return pair<int,dim_t>(1, getNumDOF());
161 case Elements:
162 return pair<int,dim_t>(ptsPerSample, getNumElements());
163 case FaceElements:
164 return pair<int,dim_t>(ptsPerSample/2, getNumFaceElements());
165 case ReducedElements:
166 return pair<int,dim_t>(1, getNumElements());
167 case ReducedFaceElements:
168 return pair<int,dim_t>(1, getNumFaceElements());
169 case Points:
170 return pair<int,dim_t>(1, m_diracPoints.size());
171 default:
172 break;
173 }
174
175 stringstream msg;
176 msg << "getDataShape: Invalid function space type " << fsType
177 << " for " << getDescription();
178 throw ValueError(msg.str());
179 }
180
181 string RipleyDomain::showTagNames() const
182 {
183 stringstream ret;
184 TagMap::const_iterator it;
185 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
186 if (it!=m_tagMap.begin()) ret << ", ";
187 ret << it->first;
188 }
189 return ret.str();
190 }
191
192 bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
193 {
194 /*
195 The idea is to use equivalence classes (i.e. types which can be
196 interpolated back and forth):
197 class 0: DOF <-> Nodes
198 class 1: ReducedDOF <-> ReducedNodes
199 class 2: Points
200 class 3: Elements
201 class 4: ReducedElements
202 class 5: FaceElements
203 class 6: ReducedFaceElements
204
205 There is also a set of lines. Interpolation is possible down a line but not
206 between lines.
207 class 0 and 1 belong to all lines so aren't considered.
208 line 0: class 2
209 line 1: class 3,4
210 line 2: class 5,6
211
212 For classes with multiple members (eg class 1) we have vars to record if
213 there is at least one instance. e.g. hasnodes is true if we have at least
214 one instance of Nodes.
215 */
216 if (fs.empty())
217 return false;
218 vector<bool> hasclass(7, false);
219 vector<int> hasline(3, 0);
220 bool hasnodes=false;
221 bool hasrednodes=false;
222 for (size_t i=0; i<fs.size(); ++i) {
223 switch (fs[i]) {
224 case Nodes: hasnodes=true; // fall through
225 case DegreesOfFreedom:
226 hasclass[0]=true;
227 break;
228 case ReducedNodes: hasrednodes=true; // fall through
229 case ReducedDegreesOfFreedom:
230 hasclass[1]=true;
231 break;
232 case Points:
233 hasline[0]=1;
234 hasclass[2]=true;
235 break;
236 case Elements:
237 hasclass[3]=true;
238 hasline[1]=1;
239 break;
240 case ReducedElements:
241 hasclass[4]=true;
242 hasline[1]=1;
243 break;
244 case FaceElements:
245 hasclass[5]=true;
246 hasline[2]=1;
247 break;
248 case ReducedFaceElements:
249 hasclass[6]=true;
250 hasline[2]=1;
251 break;
252 default:
253 return false;
254 }
255 }
256 int numLines=hasline[0]+hasline[1]+hasline[2];
257
258 // fail if we have more than one leaf group
259 // = there are at least two branches we can't interpolate between
260 if (numLines > 1)
261 return false;
262 else if (numLines==1) {
263 // we have points
264 if (hasline[0]==1)
265 resultcode=Points;
266 else if (hasline[1]==1) {
267 if (hasclass[4])
268 resultcode=ReducedElements;
269 else
270 resultcode=Elements;
271 } else { // hasline[2]==1
272 if (hasclass[6])
273 resultcode=ReducedFaceElements;
274 else
275 resultcode=FaceElements;
276 }
277 } else { // numLines==0
278 if (hasclass[1])
279 // something from class 1
280 resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
281 else // something from class 0
282 resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
283 }
284 return true;
285 }
286
287 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
288 int fsType_target) const
289 {
290 if (!isValidFunctionSpaceType(fsType_target)) {
291 stringstream msg;
292 msg << "probeInterpolationOnDomain: Invalid function space type "
293 << fsType_target << " for " << getDescription();
294 throw ValueError(msg.str());
295 }
296
297 switch (fsType_source) {
298 case Nodes:
299 case DegreesOfFreedom:
300 return true;
301 case ReducedNodes:
302 case ReducedDegreesOfFreedom:
303 return (fsType_target != Nodes &&
304 fsType_target != DegreesOfFreedom);
305 case Elements:
306 case ReducedElements:
307 return (fsType_target==Elements ||
308 fsType_target==ReducedElements);
309 case FaceElements:
310 case ReducedFaceElements:
311 return (fsType_target==FaceElements ||
312 fsType_target==ReducedFaceElements);
313 case Points:
314 return (fsType_target==fsType_source);
315
316 default: {
317 stringstream msg;
318 msg << "probeInterpolationOnDomain: Invalid function space type "
319 << fsType_source << " for " << getDescription();
320 throw ValueError(msg.str());
321 }
322 }
323 }
324
325 signed char RipleyDomain::preferredInterpolationOnDomain(int fsType_source,
326 int fsType_target) const
327 {
328 if (!isValidFunctionSpaceType(fsType_target)) {
329 stringstream msg;
330 msg << "preferredInterpolationOnDomain: Invalid function space type "
331 << fsType_target << " for " << getDescription();
332 throw ValueError(msg.str());
333 }
334
335 if (fsType_source==fsType_target) {
336 return 1;
337 }
338 // There is a complication here in that Nodes and DegreesOfFreedom
339 // can be interpolated to anything, so we need to handle the
340 // reverse case for them specially
341
342 if ((fsType_target==Nodes) || (fsType_target==DegreesOfFreedom)) {
343 return -1; // reverse interpolation
344 }
345
346 switch (fsType_source) {
347 case Nodes:
348 case DegreesOfFreedom:
349 return 1;
350 case ReducedNodes:
351 case ReducedDegreesOfFreedom:
352 return (fsType_target != Nodes &&
353 fsType_target != DegreesOfFreedom) ? -1 : 0;
354 case Elements:
355 return (fsType_target==ReducedElements) ? 1 : 0;
356 case ReducedElements:
357 return (fsType_target==Elements) ? -1 : 0;
358 case FaceElements:
359 return (fsType_target==ReducedFaceElements) ? 1 : 0;
360 case ReducedFaceElements:
361 return (fsType_target==FaceElements) ? -1 : 0;
362 case Points:
363 return false; // other case caught by the if above
364 default: {
365 stringstream msg;
366 msg << "probeInterpolationOnDomain: Invalid function space type "
367 << fsType_source << " for " << getDescription();
368 throw ValueError(msg.str());
369 }
370 }
371 }
372
373 void RipleyDomain::interpolateOnDomain(escript::Data& target,
374 const escript::Data& in) const
375 {
376 const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
377 const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
378 if (inDomain != *this)
379 throw ValueError("Illegal domain of interpolant");
380 if (targetDomain != *this)
381 throw ValueError("Illegal domain of interpolation target");
382 if (target.isComplex() != in.isComplex())
383 throw ValueError("Complexity of input and output must match");
384
385 stringstream msg;
386 msg << "interpolateOnDomain() not implemented for function space "
387 << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
388 << " -> "
389 << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
390
391 const int inFS = in.getFunctionSpace().getTypeCode();
392 const int outFS = target.getFunctionSpace().getTypeCode();
393
394 // simplest case: 1:1 copy
395 if (inFS==outFS) {
396 if (in.isComplex())
397 copyData<cplx_t>(target, in);
398 else
399 copyData<real_t>(target, in);
400 // not allowed: reduced nodes/DOF->non-reduced nodes/DOF
401 } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
402 && (outFS==Nodes || outFS==DegreesOfFreedom)) {
403 throw ValueError("interpolateOnDomain: Cannot interpolate "
404 "reduced data to non-reduced data.");
405 } else if ((inFS==Elements && outFS==ReducedElements)
406 || (inFS==FaceElements && outFS==ReducedFaceElements)) {
407 if (in.actsExpanded()) {
408 if (in.isComplex())
409 averageData<cplx_t>(target, in);
410 else
411 averageData<real_t>(target, in);
412 } else {
413 if (in.isComplex())
414 copyData<cplx_t>(target, in);
415 else
416 copyData<real_t>(target, in);
417 }
418 } else if ((inFS==ReducedElements && outFS==Elements)
419 || (inFS==ReducedFaceElements && outFS==FaceElements)) {
420 if (in.isComplex())
421 multiplyData<cplx_t>(target, in);
422 else
423 multiplyData<real_t>(target, in);
424 } else {
425 switch (inFS) {
426 case Nodes:
427 case ReducedNodes:
428 switch (outFS) {
429 case DegreesOfFreedom:
430 case ReducedDegreesOfFreedom:
431 if (getMPISize()==1) {
432 if (in.isComplex())
433 copyData<cplx_t>(target, in);
434 else
435 copyData<real_t>(target, in);
436 } else {
437 if (in.isComplex())
438 throw NotImplementedError("nodesToDOF not implemented for complex Data");
439 else
440 nodesToDOF(target, in);
441 }
442 break;
443
444 case Nodes:
445 case ReducedNodes:
446 if (in.isComplex())
447 copyData<cplx_t>(target, in);
448 else
449 copyData<real_t>(target, in);
450 break;
451 case Elements:
452 interpolateNodesOnElements(target, in, false);
453 break;
454
455 case ReducedElements:
456 interpolateNodesOnElements(target, in, true);
457 break;
458
459 case FaceElements:
460 interpolateNodesOnFaces(target, in, false);
461 break;
462
463 case ReducedFaceElements:
464 interpolateNodesOnFaces(target, in, true);
465 break;
466 case Points:
467 {
468 const dim_t numComp = in.getDataPointSize();
469 const int nDirac = m_diracPoints.size();
470 target.requireWrite();
471 #pragma omp parallel for
472 for (int i = 0; i < nDirac; i++) {
473 const double* src = in.getSampleDataRO(m_diracPoints[i].node);
474 copy(src, src+numComp, target.getSampleDataRW(i));
475 }
476 }
477 break;
478 default:
479 throw NotImplementedError(msg.str());
480 }
481 break;
482
483 case DegreesOfFreedom:
484 case ReducedDegreesOfFreedom:
485 switch (outFS) {
486 case Nodes:
487 case ReducedNodes:
488 if (getMPISize()==1)
489 if (in.isComplex())
490 copyData<cplx_t>(target, in);
491 else
492 copyData<real_t>(target, in);
493 else
494 if (in.isComplex())
495 dofToNodes<cplx_t>(target, in);
496 else
497 dofToNodes<real_t>(target, in);
498 break;
499
500 case DegreesOfFreedom:
501 case ReducedDegreesOfFreedom:
502 if (in.isComplex())
503 copyData<cplx_t>(target, in);
504 else
505 copyData<real_t>(target, in);
506 break;
507
508 case Elements:
509 case ReducedElements:
510 if (getMPISize()==1) {
511 interpolateNodesOnElements(target, in, outFS==ReducedElements);
512 } else {
513 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
514 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
515 interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
516 }
517 break;
518 case FaceElements:
519 case ReducedFaceElements:
520 if (getMPISize()==1) {
521 interpolateNodesOnFaces(target, in, outFS==ReducedFaceElements);
522 } else {
523 escript::Data contIn(in, (inFS==DegreesOfFreedom ?
524 escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
525 interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
526 }
527 break;
528
529 default:
530 throw NotImplementedError(msg.str());
531 }
532 break;
533 case Points:
534 switch(outFS) {
535 case Nodes:
536 {
537 const dim_t numComp = in.getDataPointSize();
538 const int nDirac = m_diracPoints.size();
539 target.requireWrite();
540 #pragma omp parallel for
541 for (int i = 0; i < nDirac; i++) {
542 const double* src = in.getSampleDataRO(i);
543 copy(src, src+numComp, target.getSampleDataRW(m_diracPoints[i].node));
544 }
545 }
546
547 }
548 break;
549 default:
550 throw NotImplementedError(msg.str());
551 }
552 }
553 }
554
555 escript::Data RipleyDomain::getX() const
556 {
557 return escript::continuousFunction(*this).getX();
558 }
559
560 escript::Data RipleyDomain::getNormal() const
561 {
562 return escript::functionOnBoundary(*this).getNormal();
563 }
564
565 escript::Data RipleyDomain::getSize() const
566 {
567 return escript::function(*this).getSize();
568 }
569
570 void RipleyDomain::setToX(escript::Data& arg) const
571 {
572 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
573 *(arg.getFunctionSpace().getDomain()));
574 if (argDomain != *this)
575 throw ValueError("setToX: Illegal domain of data point locations");
576 if (!arg.isExpanded())
577 throw ValueError("setToX: Expanded Data object expected");
578
579 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
580 assembleCoordinates(arg);
581 } else {
582 // interpolate the result
583 escript::Data contData=escript::Vector(0., escript::continuousFunction(*this), true);
584 assembleCoordinates(contData);
585 interpolateOnDomain(arg, contData);
586 }
587 }
588
589 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
590 {
591 const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
592 *(arg.getFunctionSpace().getDomain()));
593 if (argDomain != *this)
594 throw ValueError("setToGradient: Illegal domain of gradient argument");
595 const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
596 *(grad.getFunctionSpace().getDomain()));
597 if (gradDomain != *this)
598 throw ValueError("setToGradient: Illegal domain of gradient");
599
600 switch (grad.getFunctionSpace().getTypeCode()) {
601 case Elements:
602 case ReducedElements:
603 case FaceElements:
604 case ReducedFaceElements:
605 break;
606 default: {
607 stringstream msg;
608 msg << "setToGradient: not supported for "
609 << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
610 throw ValueError(msg.str());
611 }
612 }
613
614 switch (arg.getFunctionSpace().getTypeCode()) {
615 case DegreesOfFreedom:
616 case ReducedDegreesOfFreedom:
617 case Nodes:
618 case ReducedNodes:
619 break;
620 default: {
621 throw ValueError("setToGradient: only supported for nodal input data");
622 }
623 }
624
625 if (getMPISize() > 1) {
626 if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
627 escript::Data contArg(arg, escript::continuousFunction(*this));
628 assembleGradient(grad, contArg);
629 } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
630 escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
631 assembleGradient(grad, contArg);
632 } else {
633 assembleGradient(grad, arg);
634 }
635 } else {
636 assembleGradient(grad, arg);
637 }
638 }
639
640 void RipleyDomain::setToIntegrals(vector<real_t>& integrals, const escript::Data& arg) const
641 {
642 setToIntegralsWorker<real_t>(integrals, arg);
643 }
644
645 void RipleyDomain::setToIntegrals(vector<cplx_t>& integrals, const escript::Data& arg) const
646 {
647 setToIntegralsWorker<cplx_t>(integrals, arg);
648 }
649
650 template<typename Scalar>
651 void RipleyDomain::setToIntegralsWorker(std::vector<Scalar>& integrals,
652 const escript::Data& arg) const
653 {
654 const RipleyDomain& argDomain = dynamic_cast<const RipleyDomain&>(
655 *(arg.getFunctionSpace().getDomain()));
656 if (argDomain != *this)
657 throw ValueError("setToIntegrals: illegal domain of integration kernel");
658
659 switch (arg.getFunctionSpace().getTypeCode()) {
660 case Nodes:
661 case ReducedNodes:
662 case DegreesOfFreedom:
663 case ReducedDegreesOfFreedom:
664 {
665 escript::Data funcArg(arg, escript::function(*this));
666 assembleIntegrate(integrals, funcArg);
667 }
668 break;
669 case Elements:
670 case ReducedElements:
671 case FaceElements:
672 case ReducedFaceElements:
673 assembleIntegrate(integrals, arg);
674 break;
675 default: {
676 stringstream msg;
677 msg << "setToIntegrals: not supported for "
678 << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
679 throw ValueError(msg.str());
680 }
681 }
682 }
683
684 bool RipleyDomain::isCellOriented(int fsType) const
685 {
686 switch(fsType) {
687 case Nodes:
688 case ReducedNodes:
689 case DegreesOfFreedom:
690 case ReducedDegreesOfFreedom:
691 return false;
692 case Elements:
693 case FaceElements:
694 case Points:
695 case ReducedElements:
696 case ReducedFaceElements:
697 return true;
698 default:
699 break;
700 }
701 stringstream msg;
702 msg << "isCellOriented: invalid function space type " << fsType
703 << " on " << getDescription();
704 throw ValueError(msg.str());
705 }
706
707 bool RipleyDomain::canTag(int fsType) const
708 {
709 switch(fsType) {
710 case Nodes:
711 case Elements:
712 case ReducedElements:
713 case FaceElements:
714 case Points:
715 case ReducedFaceElements:
716 return true;
717 case DegreesOfFreedom:
718 case ReducedDegreesOfFreedom:
719 case ReducedNodes:
720 return false;
721 default:
722 break;
723 }
724 stringstream msg;
725 msg << "canTag: invalid function space type " << fsType << " on "
726 << getDescription();
727 throw ValueError(msg.str());
728 }
729
730 void RipleyDomain::setTags(int fsType, int newTag, const escript::Data& mask) const
731 {
732 vector<int>* target=NULL;
733 dim_t num=0;
734
735 switch(fsType) {
736 case Nodes:
737 num=getNumNodes();
738 target=&m_nodeTags;
739 break;
740 case Elements:
741 case ReducedElements:
742 num=getNumElements();
743 target=&m_elementTags;
744 break;
745 case FaceElements:
746 case ReducedFaceElements:
747 num=getNumFaceElements();
748 target=&m_faceTags;
749 break;
750 default: {
751 stringstream msg;
752 msg << "setTags: invalid function space type " << fsType;
753 throw ValueError(msg.str());
754 }
755 }
756 if (target->size() != num) {
757 target->assign(num, -1);
758 }
759
760 #pragma omp parallel for
761 for (index_t i=0; i<num; i++) {
762 if (mask.getSampleDataRO(i)[0] > 0) {
763 (*target)[i]=newTag;
764 }
765 }
766 updateTagsInUse(fsType);
767 }
768
769 int RipleyDomain::getTagFromSampleNo(int fsType, dim_t sampleNo) const
770 {
771 switch(fsType) {
772 case Nodes:
773 if (m_nodeTags.size() > sampleNo)
774 return m_nodeTags[sampleNo];
775 break;
776 case Elements:
777 case ReducedElements:
778 if (m_elementTags.size() > sampleNo)
779 return m_elementTags[sampleNo];
780 break;
781 case Points:
782 if (m_diracPoints.size() > sampleNo)
783 return m_diracPoints[sampleNo].tag;
784 break;
785 case FaceElements:
786 case ReducedFaceElements:
787 if (m_faceTags.size() > sampleNo)
788 return m_faceTags[sampleNo];
789 break;
790 default: {
791 stringstream msg;
792 msg << "getTagFromSampleNo: invalid function space type " << fsType;
793 throw ValueError(msg.str());
794 }
795 }
796 return -1;
797 }
798
799 int RipleyDomain::getNumberOfTagsInUse(int fsType) const
800 {
801 switch(fsType) {
802 case Nodes:
803 return m_nodeTagsInUse.size();
804 case Elements:
805 case ReducedElements:
806 return m_elementTagsInUse.size();
807 case FaceElements:
808 case ReducedFaceElements:
809 return m_faceTagsInUse.size();
810 default: {
811 stringstream msg;
812 msg << "getNumberOfTagsInUse: invalid function space type "
813 << fsType;
814 throw ValueError(msg.str());
815 }
816 }
817 }
818
819 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
820 {
821 switch(fsType) {
822 case Nodes:
823 return &m_nodeTagsInUse[0];
824 case Elements:
825 case ReducedElements:
826 return &m_elementTagsInUse[0];
827 case FaceElements:
828 case ReducedFaceElements:
829 return &m_faceTagsInUse[0];
830 default: {
831 stringstream msg;
832 msg << "borrowListOfTagsInUse: invalid function space type "
833 << fsType;
834 throw ValueError(msg.str());
835 }
836 }
837 }
838
839 void RipleyDomain::Print_Mesh_Info(bool full) const
840 {
841 cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
842 << m_mpiInfo->rank << ". MPI size: " << m_mpiInfo->size << endl;
843 cout << "Number of dimensions: " << m_numDim << endl;
844 cout << "Number of elements per rank: " << getNumElements() << endl;
845 // write tags
846 if (m_tagMap.size() > 0) {
847 cout << "Tags:" << endl;
848 TagMap::const_iterator it;
849 for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
850 cout << " " << setw(5) << it->second << " " << it->first << endl;
851 }
852 }
853 }
854
855 int RipleyDomain::getSystemMatrixTypeId(const bp::object& options) const
856 {
857 const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
858 int package = sb.getPackage();
859 escript::SolverOptions method = sb.getSolverMethod();
860 #ifdef ESYS_HAVE_TRILINOS
861 bool isDirect = escript::isDirectSolver(method);
862 #endif
863
864 // the configuration of ripley should have taken care that we have either
865 // paso or trilinos so here's how we prioritize
866 #if defined(ESYS_HAVE_PASO) && defined(ESYS_HAVE_TRILINOS)
867 // we have Paso & Trilinos so use Trilinos for parallel direct solvers
868 if (package == escript::SO_DEFAULT) {
869 if ((method == escript::SO_METHOD_DIRECT && getMPISize() > 1)
870 || isDirect
871 || sb.isComplex()) {
872 package = escript::SO_PACKAGE_TRILINOS;
873 }
874 }
875 #endif
876 #ifdef ESYS_HAVE_PASO
877 if (package == escript::SO_DEFAULT)
878 package = escript::SO_PACKAGE_PASO;
879 #endif
880 #ifdef ESYS_HAVE_TRILINOS
881 if (package == escript::SO_DEFAULT)
882 package = escript::SO_PACKAGE_TRILINOS;
883 #endif
884 if (package == escript::SO_PACKAGE_TRILINOS) {
885 #ifdef ESYS_HAVE_TRILINOS
886 int type = (int)SMT_TRILINOS;
887 if (sb.isComplex())
888 type |= (int)SMT_COMPLEX;
889 // This is required because MueLu (AMG) and Amesos2 (direct) do not
890 // support block matrices at this point. Remove if they ever do...
891 if (sb.getPreconditioner() == escript::SO_PRECONDITIONER_AMG ||
892 sb.getPreconditioner() == escript::SO_PRECONDITIONER_ILUT ||
893 isDirect) {
894 type |= (int)SMT_UNROLL;
895 }
896 return type;
897 #else
898 throw RipleyException("Trilinos requested but not built with Trilinos.");
899 #endif
900 }
901 #ifdef ESYS_HAVE_PASO
902 if (sb.isComplex()) {
903 throw NotImplementedError("Paso does not support complex-valued matrices");
904 }
905 // in all other cases we use PASO
906 return (int)SMT_PASO | paso::SystemMatrix::getSystemMatrixTypeId(
907 method, sb.getPreconditioner(), sb.getPackage(),
908 sb.isSymmetric(), m_mpiInfo);
909 #else
910 throw RipleyException("Unable to find a working solver library!");
911 #endif
912 }
913
914 int RipleyDomain::getTransportTypeId(int solver, int preconditioner,
915 int package, bool symmetry) const
916 {
917 #ifdef ESYS_USE_PASO
918 return paso::TransportProblem::getTypeId(solver, preconditioner,
919 package, symmetry, m_mpiInfo);
920 #else
921 throw RipleyException("Transport solvers require Paso but ripley was not compiled with Paso!");
922 #endif
923 }
924
925 escript::ASM_ptr RipleyDomain::newSystemMatrix(int row_blocksize,
926 const escript::FunctionSpace& row_functionspace, int column_blocksize,
927 const escript::FunctionSpace& column_functionspace, int type) const
928 {
929 bool reduceRowOrder=false;
930 bool reduceColOrder=false;
931 // is the domain right?
932 const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
933 if (row_domain != *this)
934 throw ValueError("newSystemMatrix: domain of row function space does not match the domain of matrix generator");
935 const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
936 if (col_domain != *this)
937 throw ValueError("newSystemMatrix: domain of column function space does not match the domain of matrix generator");
938 // is the function space type right?
939 if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
940 reduceRowOrder=true;
941 else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
942 throw ValueError("newSystemMatrix: illegal function space type for system matrix rows");
943 if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
944 reduceColOrder=true;
945 else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
946 throw ValueError("newSystemMatrix: illegal function space type for system matrix columns");
947 // are block sizes identical?
948 if (row_blocksize != column_blocksize)
949 throw ValueError("newSystemMatrix: row/column block sizes must be equal");
950 // are function spaces equal
951 if (reduceRowOrder != reduceColOrder)
952 throw ValueError("newSystemMatrix: row/column function spaces must be equal");
953
954 // generate matrix
955 //if (reduceRowOrder || reduceColOrder)
956 // throw NotImplementedError("newSystemMatrix: reduced order not supported");
957
958 if (type & (int)SMT_CUSP) {
959 #ifndef ESYS_HAVE_CUDA
960 throw RipleyException("eScript does not support CUDA.");
961 #endif
962 } else if (type & (int)SMT_TRILINOS) {
963 #ifdef ESYS_HAVE_TRILINOS
964 const_TrilinosGraph_ptr graph(getTrilinosGraph());
965 bool isComplex = (type & (int)SMT_COMPLEX);
966 bool unroll = (type & (int)SMT_UNROLL);
967 escript::ASM_ptr sm(new TrilinosMatrixAdapter(m_mpiInfo, row_blocksize,
968 row_functionspace, graph, isComplex, unroll));
969 return sm;
970 #else
971 throw RipleyException("newSystemMatrix: ripley was not compiled with "
972 "Trilinos support so the Trilinos solver stack cannot be used.");
973 #endif
974 } else if (type & (int)SMT_PASO) {
975 #ifdef ESYS_HAVE_PASO
976 paso::SystemMatrixPattern_ptr pattern(getPasoMatrixPattern(
977 reduceRowOrder, reduceColOrder));
978 type -= (int)SMT_PASO;
979 escript::ASM_ptr sm(new paso::SystemMatrix(type, pattern,
980 row_blocksize, column_blocksize, false, row_functionspace,
981 column_functionspace));
982 return sm;
983 #else
984 throw RipleyException("newSystemMatrix: ripley was not compiled with "
985 "Paso support so the Paso solver stack cannot be used.");
986 #endif
987 } else {
988 throw RipleyException("newSystemMatrix: unknown matrix type ID");
989 }
990 }
991
992 void RipleyDomain::addToSystem(escript::AbstractSystemMatrix& mat,
993 escript::Data& rhs, const DataMap& coefs,
994 Assembler_ptr assembler) const
995 {
996 if (isNotEmpty("d_contact", coefs) || isNotEmpty("y_contact", coefs))
997 throw ValueError(
998 "addToSystem: Ripley does not support contact elements");
999
1000 assemblePDE(&mat, rhs, coefs, assembler);
1001 assemblePDEBoundary(&mat, rhs, coefs, assembler);
1002 assemblePDEDirac(&mat, rhs, coefs, assembler);
1003 }
1004
1005 void RipleyDomain::addToSystemFromPython(escript::AbstractSystemMatrix& mat,
1006 escript::Data& rhs,
1007 const bp::list& data,
1008 Assembler_ptr assembler) const
1009 {
1010 DataMap mapping;
1011 tupleListToMap(mapping, data);
1012 addToSystem(mat, rhs, mapping, assembler);
1013 }
1014
1015 Assembler_ptr RipleyDomain::createAssemblerFromPython(const string type,
1016 const bp::list& options) const
1017 {
1018 DataMap mapping;
1019 tupleListToMap(mapping, options);
1020 return createAssembler(type, mapping);
1021 }
1022
1023 void RipleyDomain::addToRHSFromPython(escript::Data& rhs, const bp::list& data,
1024 Assembler_ptr assembler) const
1025 {
1026 DataMap mapping;
1027 tupleListToMap(mapping, data);
1028 addToRHS(rhs, mapping, assembler);
1029 }
1030
1031 void RipleyDomain::addToRHS(escript::Data& rhs, const DataMap& coefs,
1032 Assembler_ptr assembler) const
1033 {
1034 if (isNotEmpty("y_contact", coefs))
1035 throw ValueError(
1036 "addPDEToRHS: Ripley does not support contact elements");
1037
1038 if (rhs.isEmpty()) {
1039 if ((isNotEmpty("X", coefs) && isNotEmpty("du", coefs))
1040 || isNotEmpty("Y", coefs))
1041 throw ValueError(
1042 "addPDEToRHS: right hand side coefficients are provided "
1043 "but no right hand side vector given");
1044 else
1045 return;
1046 }
1047
1048 assemblePDE(NULL, rhs, coefs, assembler);
1049 assemblePDEBoundary(NULL, rhs, coefs, assembler);
1050 assemblePDEDirac(NULL, rhs, coefs, assembler);
1051 }
1052
1053 escript::ATP_ptr RipleyDomain::newTransportProblem(int blocksize,
1054 const escript::FunctionSpace& functionspace, int type) const
1055 {
1056 // is the domain right?
1057 const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));
1058 if (domain != *this)
1059 throw ValueError("newTransportProblem: domain of function space does not match the domain of transport problem generator");
1060 // is the function space type right?
1061 if (functionspace.getTypeCode() != ReducedDegreesOfFreedom &&
1062 functionspace.getTypeCode() != DegreesOfFreedom)
1063 throw ValueError("newTransportProblem: illegal function space type for transport problem");
1064
1065 #ifdef ESYS_HAVE_PASO
1066 const bool reduced = (functionspace.getTypeCode() == ReducedDegreesOfFreedom);
1067 // generate matrix
1068 paso::SystemMatrixPattern_ptr pattern(getPasoMatrixPattern(reduced, reduced));
1069 escript::ATP_ptr tp(new paso::TransportProblem(pattern, blocksize,
1070 functionspace));
1071 return tp;
1072 #else
1073 throw RipleyException("newTransportProblem: transport problems require the "
1074 "Paso library which is not available.");
1075 #endif
1076 }
1077
1078 void RipleyDomain::addPDEToTransportProblemFromPython(
1079 escript::AbstractTransportProblem& tp, escript::Data& source,
1080 const bp::list& data, Assembler_ptr assembler) const
1081 {
1082 DataMap mapping;
1083 tupleListToMap(mapping, data);
1084 addPDEToTransportProblem(tp, source, mapping, assembler);
1085 }
1086
1087 void RipleyDomain::addPDEToTransportProblem(
1088 escript::AbstractTransportProblem& tp, escript::Data& source,
1089 const DataMap& coefs, Assembler_ptr assembler) const
1090 {
1091 if (isNotEmpty("d_contact", coefs) || isNotEmpty("y_contact", coefs))
1092 throw ValueError("addPDEToTransportProblem: Ripley does not support contact elements");
1093
1094 #ifdef ESYS_HAVE_PASO
1095 paso::TransportProblem* ptp = dynamic_cast<paso::TransportProblem*>(&tp);
1096 if (!ptp)
1097 throw ValueError("addPDEToTransportProblem: Ripley only accepts Paso transport problems");
1098
1099 escript::ASM_ptr mm(ptp->borrowMassMatrix());
1100 escript::ASM_ptr tm(ptp->borrowTransportMatrix());
1101
1102 assemblePDE(mm.get(), source, coefs, assembler);
1103 assemblePDE(tm.get(), source, coefs, assembler);
1104 assemblePDEBoundary(tm.get(), source, coefs, assembler);
1105 assemblePDEDirac(tm.get(), source, coefs, assembler);
1106 #else
1107 throw RipleyException("addPDEToTransportProblem: transport problems "
1108 "require the Paso library which is not available.");
1109 #endif
1110 }
1111
1112 void RipleyDomain::addPDEToTransportProblem(
1113 escript::AbstractTransportProblem& tp, escript::Data& source,
1114 const escript::Data& M,
1115 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,
1116 const escript::Data& X,const escript::Data& Y,
1117 const escript::Data& d, const escript::Data& y,
1118 const escript::Data& d_contact,const escript::Data& y_contact,
1119 const escript::Data& d_dirac,const escript::Data& y_dirac) const
1120 {
1121 throw RipleyException("Programmer error: incorrect version of addPDEToTransportProblem called");
1122
1123 }
1124
1125 void RipleyDomain::setNewX(const escript::Data& arg)
1126 {
1127 throw NotImplementedError("setNewX(): operation not supported");
1128 }
1129
1130 //protected
1131 template<typename Scalar>
1132 void RipleyDomain::dofToNodes(escript::Data& out, const escript::Data& in) const
1133 {
1134 // expand data object if necessary
1135 const_cast<escript::Data*>(&in)->expand();
1136 const dim_t numComp = in.getDataPointSize();
1137 const dim_t numNodes = getNumNodes();
1138 const Scalar zero = static_cast<Scalar>(0);
1139 out.requireWrite();
1140 #ifdef ESYS_HAVE_TRILINOS
1141 using namespace esys_trilinos;
1142
1143 const_TrilinosGraph_ptr graph(getTrilinosGraph());
1144 Teuchos::RCP<const MapType> colMap;
1145 Teuchos::RCP<const MapType> rowMap;
1146 MapType colPointMap;
1147 MapType rowPointMap;
1148
1149 if (numComp > 1) {
1150 colPointMap = BlockVectorType<Scalar>::makePointMap(
1151 *graph->getColMap(), numComp);
1152 rowPointMap = BlockVectorType<Scalar>::makePointMap(
1153 *graph->getRowMap(), numComp);
1154 colMap = Teuchos::rcpFromRef(colPointMap);
1155 rowMap = Teuchos::rcpFromRef(rowPointMap);
1156 } else {
1157 colMap = graph->getColMap();
1158 rowMap = graph->getRowMap();
1159 }
1160
1161 const ImportType importer(rowMap, colMap);
1162 const Teuchos::ArrayView<const Scalar> localIn(
1163 in.getSampleDataRO(0, zero), in.getNumDataPoints()*numComp);
1164 Teuchos::RCP<VectorType<Scalar> > lclData = rcp(new VectorType<Scalar>(
1165 rowMap, localIn, localIn.size(), 1));
1166 Teuchos::RCP<VectorType<Scalar> > gblData = rcp(new VectorType<Scalar>(
1167 colMap, 1));
1168 gblData->doImport(*lclData, importer, Tpetra::INSERT);
1169 Teuchos::ArrayRCP<const Scalar> gblArray(gblData->getData(0));
1170
1171 #pragma omp parallel for
1172 for (index_t i = 0; i < numNodes; i++) {
1173 const Scalar* src = &gblArray[getDofOfNode(i) * numComp];
1174 copy(src, src+numComp, out.getSampleDataRW(i, zero));
1175 }
1176 #elif defined(ESYS_HAVE_PASO)
1177 paso::Coupler_ptr<Scalar> coupler(new paso::Coupler<Scalar>(m_connector,
1178 numComp, m_mpiInfo));
1179 coupler->startCollect(in.getSampleDataRO(0, zero));
1180 const dim_t numDOF = getNumDOF();
1181 const Scalar* buffer = coupler->finishCollect();
1182
1183 #pragma omp parallel for
1184 for (index_t i = 0; i < numNodes; i++) {
1185 const index_t dof = getDofOfNode(i);
1186 const Scalar* src = (dof < numDOF ? in.getSampleDataRO(dof, zero)
1187 : &buffer[(dof - numDOF) * numComp]);
1188 copy(src, src+numComp, out.getSampleDataRW(i, zero));
1189 }
1190 #endif // ESYS_HAVE_PASO
1191 }
1192
1193 //protected
1194 template<typename Scalar>
1195 void RipleyDomain::copyData(escript::Data& out, const escript::Data& in) const
1196 {
1197 const dim_t numComp = in.getDataPointSize();
1198 const dim_t numSamples = in.getNumSamples();
1199 const Scalar zero = static_cast<Scalar>(0);
1200 out.requireWrite();
1201 #pragma omp parallel for
1202 for (index_t i = 0; i < numSamples; i++) {
1203 const Scalar* src = in.getSampleDataRO(i, zero);
1204 copy(src, src+numComp, out.getSampleDataRW(i, zero));
1205 }
1206 }
1207
1208 //protected
1209 template<typename Scalar>
1210 void RipleyDomain::averageData(escript::Data& out, const escript::Data& in) const
1211 {
1212 const dim_t numComp = in.getDataPointSize();
1213 const dim_t dpp = in.getNumDataPointsPerSample();
1214 const dim_t numSamples = in.getNumSamples();
1215 const Scalar zero = static_cast<Scalar>(0);
1216 out.requireWrite();
1217 #pragma omp parallel for
1218 for (index_t i = 0; i < numSamples; i++) {
1219 const Scalar* src = in.getSampleDataRO(i, zero);
1220 Scalar* dest = out.getSampleDataRW(i, zero);
1221 for (index_t c = 0; c < numComp; c++) {
1222 Scalar res = zero;
1223 for (index_t q = 0; q < dpp; q++)
1224 res += src[c+q*numComp];
1225 *dest++ = res / static_cast<real_t>(dpp);
1226 }
1227 }
1228 }
1229
1230 //protected
1231 template<typename Scalar>
1232 void RipleyDomain::multiplyData(escript::Data& out, const escript::Data& in) const
1233 {
1234 const dim_t numComp = in.getDataPointSize();
1235 const dim_t dpp = out.getNumDataPointsPerSample();
1236 const dim_t numSamples = in.getNumSamples();
1237 const Scalar zero = static_cast<Scalar>(0);
1238 out.requireWrite();
1239 #pragma omp parallel for
1240 for (index_t i = 0; i < numSamples; i++) {
1241 const Scalar* src = in.getSampleDataRO(i, zero);
1242 Scalar* dest = out.getSampleDataRW(i, zero);
1243 for (index_t c = 0; c < numComp; c++) {
1244 for (index_t q = 0; q < dpp; q++)
1245 dest[c+q*numComp] = src[c];
1246 }
1247 }
1248 }
1249
1250 //protected
1251 void RipleyDomain::updateTagsInUse(int fsType) const
1252 {
1253 vector<int>* tagsInUse=NULL;
1254 const vector<int>* tags=NULL;
1255 switch(fsType) {
1256 case Nodes:
1257 tags=&m_nodeTags;
1258 tagsInUse=&m_nodeTagsInUse;
1259 break;
1260 case Elements:
1261 case ReducedElements:
1262 tags=&m_elementTags;
1263 tagsInUse=&m_elementTagsInUse;
1264 break;
1265 case FaceElements:
1266 case ReducedFaceElements:
1267 tags=&m_faceTags;
1268 tagsInUse=&m_faceTagsInUse;
1269 break;
1270 case Points:
1271 throw NotImplementedError("updateTagsInUse for Ripley dirac points"
1272 " not supported");
1273 default:
1274 return;
1275 }
1276
1277 // gather global unique tag values from tags into tagsInUse
1278 tagsInUse->clear();
1279 int lastFoundValue = numeric_limits<int>::min();
1280 int minFoundValue, local_minFoundValue;
1281 const int numTags = tags->size();
1282
1283 while (true) {
1284 // find smallest value bigger than lastFoundValue
1285 minFoundValue = numeric_limits<int>::max();
1286 #pragma omp parallel private(local_minFoundValue)
1287 {
1288 local_minFoundValue = minFoundValue;
1289 #pragma omp for schedule(static) nowait
1290 for (int i = 0; i < numTags; i++) {
1291 const int v = (*tags)[i];
1292 if ((v > lastFoundValue) && (v < local_minFoundValue))
1293 local_minFoundValue = v;
1294 }
1295 #pragma omp critical
1296 {
1297 if (local_minFoundValue < minFoundValue)
1298 minFoundValue = local_minFoundValue;
1299 }
1300 }
1301 #ifdef ESYS_MPI
1302 local_minFoundValue = minFoundValue;
1303 MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
1304 #endif
1305
1306 // if we found a new value add it to the tagsInUse vector
1307 if (minFoundValue < numeric_limits<int>::max()) {
1308 tagsInUse->push_back(minFoundValue);
1309 lastFoundValue = minFoundValue;
1310 } else
1311 break;
1312 }
1313 }
1314
1315 #ifdef ESYS_HAVE_PASO
1316 void RipleyDomain::createPasoConnector(const RankVector& neighbour,
1317 const IndexVector& offsetInSharedSend,
1318 const IndexVector& offsetInSharedRecv,
1319 const IndexVector& sendShared,
1320 const IndexVector& recvShared)
1321 {
1322 const index_t* sendPtr = neighbour.empty() ? NULL : &sendShared[0];
1323 const index_t* recvPtr = neighbour.empty() ? NULL : &recvShared[0];
1324 paso::SharedComponents_ptr snd_shcomp(new paso::SharedComponents(
1325 getNumDOF(), neighbour, sendPtr, offsetInSharedSend));
1326 paso::SharedComponents_ptr rcv_shcomp(new paso::SharedComponents(
1327 getNumDOF(), neighbour, recvPtr, offsetInSharedRecv));
1328 m_connector.reset(new paso::Connector(snd_shcomp, rcv_shcomp));
1329 }
1330
1331 //protected
1332 paso::Pattern_ptr RipleyDomain::createPasoPattern(
1333 const vector<IndexVector>& indices, dim_t N) const
1334 {
1335 // paso will manage the memory
1336 const dim_t M = indices.size();
1337 index_t* ptr = new index_t[M+1];
1338 ptr[0] = 0;
1339 for (index_t i=0; i < M; i++) {
1340 ptr[i+1] = ptr[i]+indices[i].size();
1341 }
1342
1343 index_t* index = new index_t[ptr[M]];
1344
1345 #pragma omp parallel for
1346 for (index_t i=0; i < M; i++) {
1347 copy(indices[i].begin(), indices[i].end(), &index[ptr[i]]);
1348 }
1349
1350 return paso::Pattern_ptr(new paso::Pattern(MATRIX_FORMAT_DEFAULT, M, N, ptr, index));
1351 }
1352 #endif // ESYS_HAVE_PASO
1353
1354 #ifdef ESYS_HAVE_TRILINOS
1355 //protected
1356 esys_trilinos::const_TrilinosGraph_ptr RipleyDomain::createTrilinosGraph(
1357 const IndexVector& myRows,
1358 const IndexVector& myColumns) const
1359 {
1360 using namespace esys_trilinos;
1361
1362 const dim_t numMatrixRows = getNumDOF();
1363
1364 TrilinosMap_ptr rowMap(new MapType(getNumDataPointsGlobal(), myRows,
1365 0, TeuchosCommFromEsysComm(m_mpiInfo->comm)));
1366
1367 IndexVector columns(getNumNodes());
1368 // order is important - our columns (=myRows) come first, followed by
1369 // shared ones (=node Id for non-DOF)
1370 #pragma omp parallel for
1371 for (size_t i=0; i<columns.size(); i++) {
1372 columns[getDofOfNode(i)] = myColumns[i];
1373 }
1374 TrilinosMap_ptr colMap(new MapType(getNumDataPointsGlobal(), columns,
1375 0, TeuchosCommFromEsysComm(m_mpiInfo->comm)));
1376
1377 // now build CSR arrays (rowPtr and colInd)
1378 const vector<IndexVector>& conns(getConnections(true));
1379 Teuchos::ArrayRCP<size_t> rowPtr(numMatrixRows+1);
1380 for (size_t i=0; i < numMatrixRows; i++) {
1381 rowPtr[i+1] = rowPtr[i] + conns[i].size();
1382 }
1383
1384 Teuchos::ArrayRCP<LO> colInd(rowPtr[numMatrixRows]);
1385
1386 #pragma omp parallel for
1387 for (index_t i=0; i < numMatrixRows; i++) {
1388 copy(conns[i].begin(), conns[i].end(), &colInd[rowPtr[i]]);
1389 }
1390
1391 TrilinosGraph_ptr graph(new GraphType(rowMap, colMap, rowPtr, colInd));
1392 Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::parameterList();
1393 params->set("Optimize Storage", true);
1394 graph->fillComplete(rowMap, rowMap, params);
1395 return graph;
1396 }
1397 #endif // ESYS_HAVE_TRILINOS
1398
1399 //protected
1400 template<>
1401 void RipleyDomain::addToSystemMatrix<real_t>(escript::AbstractSystemMatrix* mat,
1402 const IndexVector& nodes, dim_t numEq,
1403 const DoubleVector& array) const
1404 {
1405 #ifdef ESYS_HAVE_PASO
1406 paso::SystemMatrix* psm = dynamic_cast<paso::SystemMatrix*>(mat);
1407 if (psm) {
1408 addToPasoMatrix(psm, nodes, numEq, array);
1409 return;
1410 }
1411 #endif
1412 #ifdef ESYS_HAVE_CUDA
1413 SystemMatrix* rsm = dynamic_cast<SystemMatrix*>(mat);
1414 if (rsm) {
1415 rsm->add(nodes, array);
1416 return;
1417 }
1418 #endif
1419 #ifdef ESYS_HAVE_TRILINOS
1420 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1421 if (tm) {
1422 tm->add(nodes, array);
1423 return;
1424 }
1425 #endif
1426 throw RipleyException("addToSystemMatrix: unknown system matrix type");
1427 }
1428
1429 //protected
1430 template<>
1431 void RipleyDomain::addToSystemMatrix<cplx_t>(escript::AbstractSystemMatrix* mat,
1432 const IndexVector& nodes, dim_t numEq,
1433 const vector<cplx_t>& array) const
1434 {
1435 #ifdef ESYS_HAVE_TRILINOS
1436 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1437 if (tm) {
1438 tm->add(nodes, array);
1439 return;
1440 }
1441 #endif
1442 throw RipleyException("addToSystemMatrix: only Trilinos matrices support "
1443 "complex-valued assembly!");
1444 }
1445
1446 #ifdef ESYS_HAVE_PASO
1447 //private
1448 void RipleyDomain::addToPasoMatrix(paso::SystemMatrix* mat,
1449 const IndexVector& nodes, dim_t numEq,
1450 const vector<double>& array) const
1451 {
1452 const dim_t numMyCols = mat->pattern->mainPattern->numInput;
1453 const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
1454 const dim_t numSubblocks_Eq = numEq / mat->row_block_size;
1455 const dim_t numSubblocks_Sol = numEq / mat->col_block_size;
1456
1457 const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
1458 const index_t* mainBlock_index = mat->mainBlock->pattern->index;
1459 double* mainBlock_val = mat->mainBlock->val;
1460 const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
1461 const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
1462 double* col_coupleBlock_val = mat->col_coupleBlock->val;
1463 const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
1464 const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
1465 double* row_coupleBlock_val = mat->row_coupleBlock->val;
1466 index_t offset=(mat->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
1467
1468 #define UPDATE_BLOCK(VAL) do {\
1469 for (dim_t ic=0; ic<mat->col_block_size; ++ic) {\
1470 const dim_t i_Sol=ic+mat->col_block_size*l_col;\
1471 for (dim_t ir=0; ir<mat->row_block_size; ++ir) {\
1472 const dim_t i_Eq=ir+mat->row_block_size*l_row;\
1473 VAL[k*mat->block_size+ir+mat->row_block_size*ic]\
1474 += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, numEq, numEq, nodes.size())];\
1475 }\
1476 }\
1477 } while(0)
1478
1479 if (mat->type & MATRIX_FORMAT_CSC) {
1480 for (dim_t k_Sol = 0; k_Sol < nodes.size(); ++k_Sol) {
1481 // down columns of array
1482 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1483 const dim_t i_col = nodes[k_Sol]*numSubblocks_Sol+l_col;
1484 if (i_col < numMyCols) {
1485 for (dim_t k_Eq = 0; k_Eq < nodes.size(); ++k_Eq) {
1486 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1487 const dim_t i_row = nodes[k_Eq]*numSubblocks_Eq+l_row+offset;
1488 if (i_row < numMyRows+offset) {
1489 for (dim_t k = mainBlock_ptr[i_col]-offset; k < mainBlock_ptr[i_col+1]-offset; ++k) {
1490 if (mainBlock_index[k] == i_row) {
1491 UPDATE_BLOCK(mainBlock_val);
1492 break;
1493 }
1494 }
1495 } else {
1496 for (dim_t k = col_coupleBlock_ptr[i_col]-offset; k < col_coupleBlock_ptr[i_col+1]-offset; ++k) {
1497 if (row_coupleBlock_index[k] == i_row - numMyRows) {
1498 UPDATE_BLOCK(row_coupleBlock_val);
1499 break;
1500 }
1501 }
1502 }
1503 }
1504 }
1505 } else {
1506 for (dim_t k_Eq = 0; k_Eq < nodes.size(); ++k_Eq) {
1507 // across rows of array
1508 for (dim_t l_row=0; l_row<numSubblocks_Eq; ++l_row) {
1509 const dim_t i_row = nodes[k_Eq]*numSubblocks_Eq+l_row+offset;
1510 if (i_row < numMyRows+offset) {
1511 for (dim_t k = col_coupleBlock_ptr[i_col-numMyCols]-offset;
1512 k < col_coupleBlock_ptr[i_col-numMyCols+1]-offset; ++k)
1513 {
1514 if (col_coupleBlock_index[k] == i_row) {
1515 UPDATE_BLOCK(col_coupleBlock_val);
1516 break;
1517 }
1518 }
1519 }
1520 }
1521 }
1522 }
1523 }
1524 }
1525 } else {
1526 for (dim_t k_Eq = 0; k_Eq < nodes.size(); ++k_Eq) {
1527 // down columns of array
1528 for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
1529 const dim_t i_row = nodes[k_Eq]*numSubblocks_Eq+l_row;
1530 // only look at the matrix rows stored on this processor
1531 if (i_row < numMyRows) {
1532 for (dim_t k_Sol = 0; k_Sol < nodes.size(); ++k_Sol) {
1533 for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
1534 const dim_t i_col = nodes[k_Sol]*numSubblocks_Sol+l_col+offset;
1535 if (i_col < numMyCols+offset) {
1536 for (dim_t k = mainBlock_ptr[i_row]-offset; k < mainBlock_ptr[i_row+1]-offset; ++k) {
1537 if (mainBlock_index[k] == i_col) {
1538 UPDATE_BLOCK(mainBlock_val);
1539 break;
1540 }
1541 }
1542 } else {
1543 for (dim_t k = col_coupleBlock_ptr[i_row]-offset; k < col_coupleBlock_ptr[i_row+1]-offset; ++k) {
1544 if (col_coupleBlock_index[k] == i_col-numMyCols) {
1545 UPDATE_BLOCK(col_coupleBlock_val);
1546 break;
1547 }
1548 }
1549 }
1550 }
1551 }
1552 } else {
1553 for (dim_t k_Sol = 0; k_Sol < nodes.size(); ++k_Sol) {
1554 // across rows of array
1555 for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1556 const dim_t i_col = nodes[k_Sol]*numSubblocks_Sol+l_col+offset;
1557 if (i_col < numMyCols+offset) {
1558 for (dim_t k = row_coupleBlock_ptr[i_row-numMyRows]-offset;
1559 k < row_coupleBlock_ptr[i_row-numMyRows+1]-offset; ++k)
1560 {
1561 if (row_coupleBlock_index[k] == i_col) {
1562 UPDATE_BLOCK(row_coupleBlock_val);
1563 break;
1564 }
1565 }
1566 }
1567 }
1568 }
1569 }
1570 }
1571 }
1572 }
1573 #undef UPDATE_BLOCK
1574 }
1575 #endif // ESYS_HAVE_PASO
1576
1577 //private
1578 void RipleyDomain::assemblePDE(escript::AbstractSystemMatrix* mat,
1579 escript::Data& rhs, const DataMap& coefs,
1580 Assembler_ptr assembler) const
1581 {
1582 if (rhs.isEmpty() && (isNotEmpty("X", coefs) || isNotEmpty("du", coefs))
1583 && isNotEmpty("Y", coefs))
1584 throw ValueError("assemblePDE: right hand side coefficients are "
1585 "provided but no right hand side vector given");
1586
1587 vector<int> fsTypes;
1588 assembler->collateFunctionSpaceTypes(fsTypes, coefs);
1589
1590 if (fsTypes.empty()) {
1591 return;
1592 }
1593
1594 int fs=fsTypes[0];
1595 if (fs != Elements && fs != ReducedElements)
1596 throw ValueError("assemblePDE: illegal function space type for coefficients");
1597
1598 for (vector<int>::const_iterator it=fsTypes.begin()+1; it!=fsTypes.end(); it++) {
1599 if (*it != fs) {
1600 throw ValueError("assemblePDE: coefficient function spaces don't match");
1601 }
1602 }
1603
1604 int numEq, numComp;
1605 if (!mat) {
1606 if (rhs.isEmpty()) {
1607 numEq=numComp=1;
1608 } else {
1609 numEq=numComp=rhs.getDataPointSize();
1610 }
1611 } else {
1612 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
1613 throw ValueError("assemblePDE: matrix row block size and number of components of right hand side don't match");
1614 numEq = mat->getRowBlockSize();
1615 numComp = mat->getColumnBlockSize();
1616 }
1617
1618 if (numEq != numComp)
1619 throw ValueError("assemblePDE: number of equations and number of solutions don't match");
1620
1621 #ifdef ESYS_HAVE_TRILINOS
1622 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1623 if (tm) {
1624 tm->resumeFill();
1625 }
1626 #endif
1627
1628 if (numEq==1) {
1629 if (fs==ReducedElements) {
1630 assembler->assemblePDESingleReduced(mat, rhs, coefs);
1631 } else {
1632 assembler->assemblePDESingle(mat, rhs, coefs);
1633 }
1634 } else {
1635 if (fs==ReducedElements) {
1636 assembler->assemblePDESystemReduced(mat, rhs, coefs);
1637 } else {
1638 assembler->assemblePDESystem(mat, rhs, coefs);
1639 }
1640 }
1641
1642 #ifdef ESYS_HAVE_TRILINOS
1643 if (tm) {
1644 tm->fillComplete(true);
1645 }
1646 #endif
1647 }
1648
1649 //private
1650 void RipleyDomain::assemblePDEBoundary(escript::AbstractSystemMatrix* mat,
1651 escript::Data& rhs, const DataMap& coefs, Assembler_ptr assembler) const
1652 {
1653 if (rhs.isEmpty() && isNotEmpty("y", coefs))
1654 throw ValueError("assemblePDEBoundary: y provided but no right hand side vector given");
1655
1656 int fs=-1;
1657 if (isNotEmpty("d", coefs))
1658 fs=coefs.find("d")->second.getFunctionSpace().getTypeCode();
1659 if (isNotEmpty("y", coefs)) {
1660 DataMap::const_iterator iy = coefs.find("y");
1661 if (fs == -1)
1662 fs = iy->second.getFunctionSpace().getTypeCode();
1663 else if (fs != iy->second.getFunctionSpace().getTypeCode())
1664 throw ValueError("assemblePDEBoundary: coefficient function spaces don't match");
1665 }
1666 if (fs==-1) {
1667 return;
1668 }
1669
1670 if (fs != FaceElements && fs != ReducedFaceElements)
1671 throw ValueError("assemblePDEBoundary: illegal function space type for coefficients");
1672
1673 int numEq, numComp;
1674 if (!mat) {
1675 if (rhs.isEmpty()) {
1676 numEq=numComp=1;
1677 } else {
1678 numEq=numComp=rhs.getDataPointSize();
1679 }
1680 } else {
1681 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
1682 throw ValueError("assemblePDEBoundary: matrix row block size and number of components of right hand side don't match");
1683 numEq = mat->getRowBlockSize();
1684 numComp = mat->getColumnBlockSize();
1685 }
1686
1687 if (numEq != numComp)
1688 throw ValueError("assemblePDEBoundary: number of equations and number of solutions don't match");
1689
1690 #ifdef ESYS_HAVE_TRILINOS
1691 TrilinosMatrixAdapter* tm = dynamic_cast<TrilinosMatrixAdapter*>(mat);
1692 if (tm) {
1693 tm->resumeFill();
1694 }
1695 #endif
1696
1697 if (numEq==1) {
1698 if (fs==ReducedFaceElements)
1699 assembler->assemblePDEBoundarySingleReduced(mat, rhs, coefs);
1700 else
1701 assembler->assemblePDEBoundarySingle(mat, rhs, coefs);
1702 } else {
1703 if (fs==ReducedFaceElements)
1704 assembler->assemblePDEBoundarySystemReduced(mat, rhs, coefs);
1705 else
1706 assembler->assemblePDEBoundarySystem(mat, rhs, coefs);
1707 }
1708
1709 #ifdef ESYS_HAVE_TRILINOS
1710 if (tm) {
1711 tm->fillComplete(true);
1712 }
1713 #endif
1714 }
1715
1716 void RipleyDomain::assemblePDEDirac(escript::AbstractSystemMatrix* mat,
1717 escript::Data& rhs, const DataMap& coefs,
1718 Assembler_ptr assembler) const
1719 {
1720 bool yNotEmpty = isNotEmpty("y_dirac", coefs);
1721 bool dNotEmpty = isNotEmpty("d_dirac", coefs);
1722 if (!(yNotEmpty || dNotEmpty)) {
1723 return;
1724 }
1725 escript::Data d = unpackData("d_dirac", coefs);
1726 escript::Data y = unpackData("y_dirac", coefs);
1727 int nEq, nComp;
1728 if (!mat) {
1729 if (rhs.isEmpty()) {
1730 nEq=nComp=1;
1731 } else {
1732 nEq=nComp=rhs.getDataPointSize();
1733 }
1734 } else {
1735 if (!rhs.isEmpty() && rhs.getDataPointSize() != mat->getRowBlockSize())
1736 throw ValueError("assemblePDEDirac: matrix row block size "
1737 "and number of components of right hand side don't match");
1738 nEq = mat->getRowBlockSize();
1739 nComp = mat->getColumnBlockSize();
1740 }
1741
1742 rhs.requireWrite();
1743 for (int i = 0; i < m_diracPoints.size(); i++) { //only for this rank
1744 const index_t dof = getDofOfNode(m_diracPoints[i].node);
1745 if (yNotEmpty) {
1746 const double *EM_F = y.getSampleDataRO(i);
1747 double *F_p = rhs.getSampleDataRW(0);
1748 if (dof < getNumDOF()) {
1749 for (index_t eq = 0; eq < nEq; eq++) {
1750 F_p[INDEX2(eq, dof, nEq)] += EM_F[eq];
1751 }
1752 }
1753 }
1754 if (dNotEmpty) {
1755 const IndexVector rowIndex(1, dof);
1756 const double *EM_S = d.getSampleDataRO(i);
1757 vector<double> contents(EM_S, EM_S+nEq*nEq*nComp);
1758 addToSystemMatrix(mat, rowIndex, nEq, contents);
1759 }
1760 }
1761 }
1762
1763 bool RipleyDomain::probeInterpolationAcross(int fsType_source,
1764 const escript::AbstractDomain&, int fsType_target) const
1765 {
1766 return false;
1767 }
1768
1769 void RipleyDomain::interpolateAcross(escript::Data& target,
1770 const escript::Data& source) const
1771 {
1772 throw NotImplementedError("interpolateAcross() not supported");
1773 }
1774
1775 // Expecting ("gaussian", radius, sigma)
1776 bool RipleyDomain::supportsFilter(const bp::tuple& t) const
1777 {
1778 if (len(t) == 0) { // so we can handle unfiltered randoms
1779 return true;
1780 }
1781 if (len(t) != 3) {
1782 return false;
1783 }
1784 bp::extract<string> ex(t[0]);
1785 if (!ex.check() || (ex() != "gaussian")) {
1786 return false;
1787 }
1788 if (! bp::extract<unsigned int>(t[1]).check()) {
1789 return false;
1790 }
1791 return bp::extract<double>(t[2]).check();
1792 }
1793
1794 void RipleyDomain::addPoints(const vector<double>& coords,
1795 const vector<int>& tags)
1796 {
1797 for (int i = 0; i < tags.size(); i++) {
1798 dim_t node = findNode(&coords[i * m_numDim]);
1799 if (node >= 0) {
1800 m_diracPointNodeIDs.push_back(borrowSampleReferenceIDs(Nodes)[node]);
1801 DiracPoint dp;
1802 dp.node = node; //local
1803 dp.tag = tags[i];
1804 m_diracPoints.push_back(dp);
1805 }
1806 }
1807 }
1808
1809 } // end of namespace ripley
1810

  ViewVC Help
Powered by ViewVC 1.1.26