/[escript]/branches/domexper/finley/src/CPPAdapter/TransportProblemAdapter.cpp
ViewVC logotype

Contents of /branches/domexper/finley/src/CPPAdapter/TransportProblemAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3234 - (show annotations)
Mon Oct 4 01:46:30 2010 UTC (9 years, 2 months ago) by jfenwick
File size: 6253 byte(s)
Some subdirs need to have changes pulled over but all of the unit tests 
except for modellib appear to work

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26