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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3234 - (hide annotations)
Mon Oct 4 01:46:30 2010 UTC (9 years, 3 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 ksteube 1811
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jfenwick 3234 #ifdef ESYS_MPI
16 ksteube 817 #include <mpi.h>
17     #endif
18 gross 1364 #include "TransportProblemAdapter.h"
19 jgs 203 #include "SystemMatrixAdapter.h"
20 jgs 82
21     using namespace std;
22    
23     namespace finley {
24    
25 jgs 102 struct null_deleter
26     {
27     void operator()(void const *ptr) const
28     {
29     }
30     };
31    
32    
33 gross 1364 TransportProblemAdapter::TransportProblemAdapter()
34 jgs 82 {
35 gross 1364 throw FinleyAdapterException("Error - Illegal to generate default TransportProblemAdapter.");
36 jgs 82 }
37    
38 gross 1364
39 gross 2987 TransportProblemAdapter::TransportProblemAdapter(Paso_TransportProblem* transport_problem,
40     const bool useBackwardEuler,
41 gross 1364 const int block_size,
42     const escript::FunctionSpace& functionspace):
43 gross 2987 AbstractTransportProblem(useBackwardEuler, block_size, functionspace)
44 jgs 82 {
45 gross 1364 m_transport_problem.reset(transport_problem,null_deleter());
46 jgs 82 }
47    
48 gross 1364 TransportProblemAdapter::~TransportProblemAdapter()
49 jgs 82 {
50 gross 1364 if (m_transport_problem.unique()) {
51 gross 2987 Paso_TransportProblem* transp=m_transport_problem.get();
52     Paso_TransportProblem_free(transp);
53 jgs 82 }
54     }
55    
56 gross 2987 Paso_TransportProblem* TransportProblemAdapter::getPaso_TransportProblem() const
57 jgs 82 {
58 gross 1364 return m_transport_problem.get();
59 jgs 82 }
60    
61 jgs 150
62 gross 2987 void TransportProblemAdapter::setToSolution(escript::Data& out, escript::Data& u0, escript::Data& source, const double dt, boost::python::object& options) const
63 ksteube 1339 {
64 gross 2987 Paso_TransportProblem* transp=getPaso_TransportProblem();
65 jgs 150 Paso_Options paso_options;
66 gross 2474 SystemMatrixAdapter::escriptToPasoOptions(&paso_options,options);
67     options.attr("resetDiagnostics")();
68 gross 1364 if ( out.getDataPointSize() != getBlockSize()) {
69 gross 1365 throw FinleyAdapterException("solve : block size of solution does not match block size of transport problems.");
70 gross 1364 } else if ( source.getDataPointSize() != getBlockSize()) {
71 gross 1365 throw FinleyAdapterException("solve : block size of source term does not match block size of transport problems.");
72 gross 1364 } else if ( out.getFunctionSpace() != getFunctionSpace()) {
73 gross 1365 throw FinleyAdapterException("solve : function spaces of solution and of transport problem don't match.");
74 gross 1364 } else if (source.getFunctionSpace() != getFunctionSpace()) {
75 gross 1365 throw FinleyAdapterException("solve : function spaces of source term and of transport problem don't match.");
76 gross 1364 } else if (dt<=0.) {
77 gross 1365 throw FinleyAdapterException("solve : time increment dt needs to be positive.");
78 jgs 150 }
79 jgs 153 out.expand();
80 gross 1364 source.expand();
81 jfenwick 2271 out.requireWrite();
82     source.requireWrite();
83     double* out_dp=out.getSampleDataRW(0);
84 gross 2987 double* u0_dp=u0.getSampleDataRW(0);
85 jfenwick 2271 double* source_dp=source.getSampleDataRW(0);
86 gross 2987 Paso_TransportProblem_solve(transp,out_dp,dt,u0_dp,source_dp,&paso_options);
87 gross 2474 SystemMatrixAdapter::pasoToEscriptOptions(&paso_options,options);
88 jgs 150 checkPasoError();
89 jgs 82 }
90    
91 gross 1364 void TransportProblemAdapter::resetTransport() const
92 jgs 82 {
93 gross 2987 Paso_TransportProblem* transp = getPaso_TransportProblem();
94     Paso_TransportProblem_reset(transp);
95 gross 1364 checkPasoError();
96 jgs 82 }
97 gross 2197 void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r, const double factor) const
98 gross 1417 {
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 gross 2987 Paso_TransportProblem* transp=getPaso_TransportProblem();
113 jgs 82
114 gross 1417 /* 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 jfenwick 2271 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 gross 1417
129 gross 2197 if (false) {
130     cout << "v1\n";
131     Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
132     checkPasoError();
133 gross 1417
134 gross 2197 /* 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 gross 2987 Paso_TransportProblem_setUpConstraint(transp, q_dp, factor);
145 gross 2197 checkPasoError();
146 gross 2987 Paso_TransportProblem_insertConstraint(transp,r2_dp, source_dp);
147 gross 2197 checkPasoError();
148     }
149 gross 1417 }
150    
151 gross 1407 double TransportProblemAdapter::getSafeTimeStepSize() const
152     {
153 gross 2987 Paso_TransportProblem* transp=getPaso_TransportProblem();
154     double dt=Paso_TransportProblem_getSafeTimeStepSize(transp);
155 gross 1407 checkPasoError();
156     return dt;
157     }
158 woo409 757
159 gross 1859 double TransportProblemAdapter::getUnlimitedTimeStepSize() const
160     {
161     return LARGE_POSITIVE_FLOAT;
162     }
163 gross 1407
164    
165 gross 1859
166 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26