/[escript]/release/3.4.2/pasowrap/src/TransportProblemAdapter.cpp
ViewVC logotype

Contents of /release/3.4.2/pasowrap/src/TransportProblemAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4919 - (show annotations)
Wed Apr 30 06:25:55 2014 UTC (6 years, 1 month ago) by jfenwick
File size: 6453 byte(s)
Because we've got to!

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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 "TransportProblemAdapter.h"
18 #include "SystemMatrixAdapter.h"
19
20 using namespace std;
21
22 namespace paso {
23
24 TransportProblemAdapter::TransportProblemAdapter()
25 {
26 throw PasoException("Error - Illegal to generate default TransportProblemAdapter.");
27 }
28
29 TransportProblemAdapter::TransportProblemAdapter(TransportProblem_ptr tp,
30 int block_size, const escript::FunctionSpace& functionspace) :
31 AbstractTransportProblem(block_size, functionspace),
32 m_transport_problem(tp)
33 {
34 }
35
36 TransportProblem_ptr TransportProblemAdapter::getPaso_TransportProblem() const
37 {
38 return m_transport_problem;
39 }
40
41 void TransportProblemAdapter::setToSolution(escript::Data& out,
42 escript::Data& u0, escript::Data& source, double dt,
43 boost::python::object& options) const
44 {
45 Options paso_options;
46 SystemMatrixAdapter::escriptToPasoOptions(&paso_options, options);
47 options.attr("resetDiagnostics")();
48 if ( out.getDataPointSize() != getBlockSize()) {
49 throw PasoException("solve : block size of solution does not match block size of transport problems.");
50 } else if ( source.getDataPointSize() != getBlockSize()) {
51 throw PasoException("solve : block size of source term does not match block size of transport problems.");
52 } else if ( out.getFunctionSpace() != getFunctionSpace()) {
53 throw PasoException("solve : function spaces of solution and of transport problem don't match.");
54 } else if (source.getFunctionSpace() != getFunctionSpace()) {
55 throw PasoException("solve : function spaces of source term and of transport problem don't match.");
56 } else if (dt<=0.) {
57 throw PasoException("solve : time increment dt needs to be positive.");
58 }
59 out.expand();
60 source.expand();
61 u0.expand();
62 out.requireWrite();
63 source.requireWrite();
64 double* out_dp = out.getSampleDataRW(0);
65 double* u0_dp = u0.getSampleDataRW(0);
66 double* source_dp = source.getSampleDataRW(0);
67 SystemMatrixAdapter::pasoToEscriptOptions(&paso_options, options);
68 m_transport_problem->solve(out_dp, dt, u0_dp, source_dp, &paso_options);
69
70 checkPasoError();
71 }
72
73 void TransportProblemAdapter::resetTransport() const
74 {
75 m_transport_problem->reset();
76 checkPasoError();
77 }
78
79 void TransportProblemAdapter::copyConstraint(escript::Data& source,
80 escript::Data& q,
81 escript::Data& r) const
82 {
83 if (q.getDataPointSize() != getBlockSize()) {
84 throw PasoException("copyConstraint : block size does not match the number of components of constraint mask.");
85 } else if ( q.getFunctionSpace() != getFunctionSpace()) {
86 throw PasoException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
87 } else if ( r.getDataPointSize() != getBlockSize()) {
88 throw PasoException("copyConstraint : block size does not match the number of components of constraint values.");
89 } else if ( r.getFunctionSpace() != getFunctionSpace()) {
90 throw PasoException("copyConstraint : function spaces of transport problem and constraint values don't match.");
91 } else if ( source.getDataPointSize() != getBlockSize()) {
92 throw PasoException("copyConstraint : block size does not match the number of components of source.");
93 } else if ( source.getFunctionSpace() != getFunctionSpace()) {
94 throw PasoException("copyConstraint : function spaces of transport problem and source don't match.");
95 }
96
97 if (false) {
98 // r2=r where q>0, 0 elsewhere
99 escript::Data r2(0., q.getDataPointShape(), q.getFunctionSpace());
100 r2.copyWithMask(r, q);
101
102 // source -= tp->mass_matrix*r2
103 r2.expand();
104 source.expand();
105 q.expand();
106 r2.requireWrite();
107 source.requireWrite();
108 q.requireWrite();
109 double* r2_dp = r2.getSampleDataRW(0);
110 double* source_dp = source.getSampleDataRW(0);
111 double* q_dp = q.getSampleDataRW(0);
112
113 SystemMatrix_MatrixVector(-1., m_transport_problem->mass_matrix,
114 r2_dp, 1., source_dp);
115 checkPasoError();
116
117 // insert 0 rows into transport matrix
118 m_transport_problem->transport_matrix->nullifyRows(q_dp, 0.);
119 checkPasoError();
120
121 // insert 0 rows and 1 in main diagonal into mass matrix
122 m_transport_problem->mass_matrix->nullifyRowsAndCols(q_dp, q_dp, 1.);
123 checkPasoError();
124
125 source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
126 } else {
127 r.expand();
128 source.expand();
129 q.expand();
130 r.requireWrite();
131 source.requireWrite();
132 q.requireWrite();
133 double* r_dp = r.getSampleDataRW(0);
134 double* source_dp = source.getSampleDataRW(0);
135 double* q_dp = q.getSampleDataRW(0);
136 m_transport_problem->setUpConstraint(q_dp);
137 checkPasoError();
138 m_transport_problem->insertConstraint(r_dp, source_dp);
139 checkPasoError();
140 }
141 }
142
143 double TransportProblemAdapter::getSafeTimeStepSize() const
144 {
145 const double dt = m_transport_problem->getSafeTimeStepSize();
146 checkPasoError();
147 return dt;
148 }
149
150 double TransportProblemAdapter::getUnlimitedTimeStepSize() const
151 {
152 return LARGE_POSITIVE_FLOAT;
153 }
154
155 int TransportProblemAdapter::getTransportTypeId(int solver, int preconditioner,
156 int package, bool symmetry,
157 Esys_MPIInfo* mpiInfo)
158 {
159 return TransportProblem::getTypeId(
160 SystemMatrixAdapter::mapOptionToPaso(solver),
161 SystemMatrixAdapter::mapOptionToPaso(preconditioner),
162 SystemMatrixAdapter::mapOptionToPaso(package), symmetry, mpiInfo);
163 }
164
165
166 } // end of namespace
167

  ViewVC Help
Powered by ViewVC 1.1.26