/[escript]/branches/symbolic_from_3470/pasowrap/src/TransportProblemAdapter.cpp
ViewVC logotype

Annotation of /branches/symbolic_from_3470/pasowrap/src/TransportProblemAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3868 - (hide annotations)
Thu Mar 15 06:07:08 2012 UTC (7 years, 2 months ago) by caltinay
File size: 6940 byte(s)
Update to latest trunk

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

  ViewVC Help
Powered by ViewVC 1.1.26