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

Contents of /branches/symbolic_from_3470/dudley/src/CPPAdapter/TransportProblemAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3471 - (show annotations)
Tue Mar 15 04:23:54 2011 UTC (8 years, 1 month ago) by caltinay
File size: 6212 byte(s)
Branch to investigate how to incorporate sympy or other symbolic evaluation
libraries into escript.

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