/[escript]/trunk/pasowrap/src/TransportProblemAdapter.cpp
ViewVC logotype

Contents of /trunk/pasowrap/src/TransportProblemAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 6 months ago) by jfenwick
File size: 7076 byte(s)
Round 1 of copyright fixes
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 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 since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 #include "TransportProblemAdapter.h"
17 #include "SystemMatrixAdapter.h"
18
19 using namespace std;
20
21 namespace paso {
22
23 struct null_deleter
24 {
25 void operator()(void const *ptr) const
26 {
27 }
28 };
29
30 PASOWRAP_DLL_API
31 TransportProblemAdapter::TransportProblemAdapter()
32 {
33 throw PasoException("Error - Illegal to generate default TransportProblemAdapter.");
34 }
35
36
37 PASOWRAP_DLL_API
38 TransportProblemAdapter::TransportProblemAdapter(Paso_TransportProblem* transport_problem,
39 const int block_size,
40 const escript::FunctionSpace& functionspace):
41 AbstractTransportProblem(block_size, functionspace)
42 {
43 m_transport_problem.reset(transport_problem,null_deleter());
44 }
45
46 PASOWRAP_DLL_API
47 TransportProblemAdapter::~TransportProblemAdapter()
48 {
49 if (m_transport_problem.unique()) {
50 Paso_TransportProblem* transp=m_transport_problem.get();
51 Paso_TransportProblem_free(transp);
52 }
53 }
54
55 PASOWRAP_DLL_API
56 Paso_TransportProblem* TransportProblemAdapter::getPaso_TransportProblem() const
57 {
58 return m_transport_problem.get();
59 }
60
61 PASOWRAP_DLL_API
62 void TransportProblemAdapter::setToSolution(escript::Data& out, escript::Data& u0, escript::Data& source, const double dt, boost::python::object& options) const
63 {
64 Paso_TransportProblem* transp=getPaso_TransportProblem();
65 Paso_Options paso_options;
66 SystemMatrixAdapter::escriptToPasoOptions(&paso_options,options);
67 options.attr("resetDiagnostics")();
68 if ( out.getDataPointSize() != getBlockSize()) {
69 throw PasoException("solve : block size of solution does not match block size of transport problems.");
70 } else if ( source.getDataPointSize() != getBlockSize()) {
71 throw PasoException("solve : block size of source term does not match block size of transport problems.");
72 } else if ( out.getFunctionSpace() != getFunctionSpace()) {
73 throw PasoException("solve : function spaces of solution and of transport problem don't match.");
74 } else if (source.getFunctionSpace() != getFunctionSpace()) {
75 throw PasoException("solve : function spaces of source term and of transport problem don't match.");
76 } else if (dt<=0.) {
77 throw PasoException("solve : time increment dt needs to be positive.");
78 }
79 out.expand();
80 source.expand();
81 u0.expand();
82 out.requireWrite();
83 source.requireWrite();
84 double* out_dp=out.getSampleDataRW(0);
85 double* u0_dp=u0.getSampleDataRW(0);
86 double* source_dp=source.getSampleDataRW(0);
87 SystemMatrixAdapter::pasoToEscriptOptions(&paso_options,options);
88 Paso_TransportProblem_solve(transp,out_dp,dt,u0_dp,source_dp,&paso_options);
89
90 checkPasoError();
91 }
92
93 PASOWRAP_DLL_API
94 void TransportProblemAdapter::resetTransport() const
95 {
96 Paso_TransportProblem* transp = getPaso_TransportProblem();
97 Paso_TransportProblem_reset(transp);
98 checkPasoError();
99 }
100
101 PASOWRAP_DLL_API
102 void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r) const
103 {
104 if ( q.getDataPointSize() != getBlockSize()) {
105 throw PasoException("copyConstraint : block size does not match the number of components of constraint mask.");
106 } else if ( q.getFunctionSpace() != getFunctionSpace()) {
107 throw PasoException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
108 } else if ( r.getDataPointSize() != getBlockSize()) {
109 throw PasoException("copyConstraint : block size does not match the number of components of constraint values.");
110 } else if ( r.getFunctionSpace() != getFunctionSpace()) {
111 throw PasoException("copyConstraint : function spaces of transport problem and constraint values don't match.");
112 } else if ( source.getDataPointSize() != getBlockSize()) {
113 throw PasoException("copyConstraint : block size does not match the number of components of source.");
114 } else if ( source.getFunctionSpace() != getFunctionSpace()) {
115 throw PasoException("copyConstraint : function spaces of transport problem and source don't match.");
116 }
117 Paso_TransportProblem* transp=getPaso_TransportProblem();
118
119
120 if (false) {
121
122 /* r2=r where q>0, 0 elsewhere */
123 escript::Data r2(0.,q.getDataPointShape(),q.getFunctionSpace());
124 r2.copyWithMask(r,q);
125
126 /* source-=transp->mass_matrix*r2 */
127 r2.expand();
128 source.expand();
129 q.expand();
130 r2.requireWrite();
131 source.requireWrite();
132 q.requireWrite();
133 double* r2_dp=r2.getSampleDataRW(0);
134 double* source_dp=source.getSampleDataRW(0);
135 double* q_dp=q.getSampleDataRW(0);
136
137 cout << "v1\n";
138 Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
139 checkPasoError();
140
141 /* insert 0 rows into transport matrix */
142 Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
143 checkPasoError();
144
145 /* insert 0 rows amd 1 in main diagonal into mass matrix */
146 Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
147 checkPasoError();
148
149 source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
150 } else {
151 r.expand();
152 source.expand();
153 q.expand();
154 r.requireWrite();
155 source.requireWrite();
156 q.requireWrite();
157 double* r_dp=r.getSampleDataRW(0);
158 double* source_dp=source.getSampleDataRW(0);
159 double* q_dp=q.getSampleDataRW(0);
160 Paso_TransportProblem_setUpConstraint(transp, q_dp);
161 checkPasoError();
162 Paso_TransportProblem_insertConstraint(transp,r_dp, source_dp);
163 checkPasoError();
164 }
165 }
166
167 PASOWRAP_DLL_API
168 double TransportProblemAdapter::getSafeTimeStepSize() const
169 {
170 Paso_TransportProblem* transp=getPaso_TransportProblem();
171 double dt=Paso_TransportProblem_getSafeTimeStepSize(transp);
172 checkPasoError();
173 return dt;
174 }
175
176 PASOWRAP_DLL_API
177 double TransportProblemAdapter::getUnlimitedTimeStepSize() const
178 {
179 return LARGE_POSITIVE_FLOAT;
180 }
181
182 PASOWRAP_DLL_API
183 int TransportProblemAdapter::getTransportTypeId(const int solver,
184 const int preconditioner, const int package, const bool symmetry,
185 Esys_MPIInfo* mpiInfo)
186 {
187 int out=Paso_TransportProblem_getTypeId(
188 SystemMatrixAdapter::mapOptionToPaso(solver),
189 SystemMatrixAdapter::mapOptionToPaso(preconditioner),
190 SystemMatrixAdapter::mapOptionToPaso(package),
191 symmetry?1:0, mpiInfo);
192 checkPasoError();
193 return out;
194 }
195
196
197 } // end of namespace
198

  ViewVC Help
Powered by ViewVC 1.1.26