/[escript]/trunk/dudley/src/CPPAdapter/TransportProblemAdapter.cpp
ViewVC logotype

Contents of /trunk/dudley/src/CPPAdapter/TransportProblemAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3317 - (show annotations)
Thu Oct 28 00:50:41 2010 UTC (8 years, 9 months ago) by caltinay
File size: 6212 byte(s)
Removed bogus mpi.h includes, replaced others by our esysUtils wrapper
and rearranged so that the wrapper always comes before netcdf which fixes
linking problems when disabling mpi on systems with netCDF 4.x.

1
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 dudley {
20
21 struct null_deleter
22 {
23 void operator()(void const *ptr) const
24 {
25 }
26 };
27
28
29 TransportProblemAdapter::TransportProblemAdapter()
30 {
31 throw DudleyAdapterException("Error - Illegal to generate default TransportProblemAdapter.");
32 }
33
34
35 TransportProblemAdapter::TransportProblemAdapter(Paso_TransportProblem* transport_problem,
36 const bool useBackwardEuler,
37 const int block_size,
38 const escript::FunctionSpace& functionspace):
39 AbstractTransportProblem(useBackwardEuler, block_size, functionspace)
40 {
41 m_transport_problem.reset(transport_problem,null_deleter());
42 }
43
44 TransportProblemAdapter::~TransportProblemAdapter()
45 {
46 if (m_transport_problem.unique()) {
47 Paso_TransportProblem* transp=m_transport_problem.get();
48 Paso_TransportProblem_free(transp);
49 }
50 }
51
52 Paso_TransportProblem* TransportProblemAdapter::getPaso_TransportProblem() const
53 {
54 return m_transport_problem.get();
55 }
56
57
58 void TransportProblemAdapter::setToSolution(escript::Data& out, escript::Data& u0, escript::Data& source, const double dt, boost::python::object& options) const
59 {
60 Paso_TransportProblem* transp=getPaso_TransportProblem();
61 Paso_Options paso_options;
62 SystemMatrixAdapter::escriptToPasoOptions(&paso_options,options);
63 options.attr("resetDiagnostics")();
64 if ( out.getDataPointSize() != getBlockSize()) {
65 throw DudleyAdapterException("solve : block size of solution does not match block size of transport problems.");
66 } else if ( source.getDataPointSize() != getBlockSize()) {
67 throw DudleyAdapterException("solve : block size of source term does not match block size of transport problems.");
68 } else if ( out.getFunctionSpace() != getFunctionSpace()) {
69 throw DudleyAdapterException("solve : function spaces of solution and of transport problem don't match.");
70 } else if (source.getFunctionSpace() != getFunctionSpace()) {
71 throw DudleyAdapterException("solve : function spaces of source term and of transport problem don't match.");
72 } else if (dt<=0.) {
73 throw DudleyAdapterException("solve : time increment dt needs to be positive.");
74 }
75 out.expand();
76 source.expand();
77 out.requireWrite();
78 source.requireWrite();
79 double* out_dp=out.getSampleDataRW(0);
80 double* u0_dp=u0.getSampleDataRW(0);
81 double* source_dp=source.getSampleDataRW(0);
82 Paso_TransportProblem_solve(transp,out_dp,dt,u0_dp,source_dp,&paso_options);
83 SystemMatrixAdapter::pasoToEscriptOptions(&paso_options,options);
84 checkPasoError();
85 }
86
87 void TransportProblemAdapter::resetTransport() const
88 {
89 Paso_TransportProblem* transp = getPaso_TransportProblem();
90 Paso_TransportProblem_reset(transp);
91 checkPasoError();
92 }
93 void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r, const double factor) const
94 {
95 if ( q.getDataPointSize() != getBlockSize()) {
96 throw DudleyAdapterException("copyConstraint : block size does not match the number of components of constraint mask.");
97 } else if ( q.getFunctionSpace() != getFunctionSpace()) {
98 throw DudleyAdapterException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
99 } else if ( r.getDataPointSize() != getBlockSize()) {
100 throw DudleyAdapterException("copyConstraint : block size does not match the number of components of constraint values.");
101 } else if ( r.getFunctionSpace() != getFunctionSpace()) {
102 throw DudleyAdapterException("copyConstraint : function spaces of transport problem and constraint values don't match.");
103 } else if ( source.getDataPointSize() != getBlockSize()) {
104 throw DudleyAdapterException("copyConstraint : block size does not match the number of components of source.");
105 } else if ( source.getFunctionSpace() != getFunctionSpace()) {
106 throw DudleyAdapterException("copyConstraint : function spaces of transport problem and source don't match.");
107 }
108 Paso_TransportProblem* transp=getPaso_TransportProblem();
109
110 /* r2=r where q>0, 0 elsewhere */
111 escript::Data r2(0.,q.getDataPointShape(),q.getFunctionSpace());
112 r2.copyWithMask(r,q);
113
114 /* source-=transp->mass_matrix*r2 */
115 r2.expand();
116 source.expand();
117 q.expand();
118 r2.requireWrite();
119 source.requireWrite();
120 q.requireWrite();
121 double* r2_dp=r2.getSampleDataRW(0);
122 double* source_dp=source.getSampleDataRW(0);
123 double* q_dp=q.getSampleDataRW(0);
124
125 if (false) {
126 cout << "v1\n";
127 Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
128 checkPasoError();
129
130 /* insert 0 rows into transport matrix */
131 Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
132 checkPasoError();
133
134 /* insert 0 rows amd 1 in main diagonal into mass matrix */
135 Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
136 checkPasoError();
137
138 source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
139 } else {
140 Paso_TransportProblem_setUpConstraint(transp, q_dp, factor);
141 checkPasoError();
142 Paso_TransportProblem_insertConstraint(transp,r2_dp, source_dp);
143 checkPasoError();
144 }
145 }
146
147 double TransportProblemAdapter::getSafeTimeStepSize() const
148 {
149 Paso_TransportProblem* transp=getPaso_TransportProblem();
150 double dt=Paso_TransportProblem_getSafeTimeStepSize(transp);
151 checkPasoError();
152 return dt;
153 }
154
155 double TransportProblemAdapter::getUnlimitedTimeStepSize() const
156 {
157 return LARGE_POSITIVE_FLOAT;
158 }
159
160
161
162 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26