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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1859 - (show annotations)
Wed Oct 8 03:03:37 2008 UTC (10 years, 7 months ago) by gross
File size: 6310 byte(s)
first version of testing for transport solver.
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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 PASO_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_FCTransportProblem* transport_problem,
40 const double theta,
41 const int block_size,
42 const escript::FunctionSpace& functionspace):
43 AbstractTransportProblem(theta, 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_FCTransportProblem* transp=m_transport_problem.get();
52 Paso_FCTransportProblem_free(transp);
53 }
54 }
55
56 Paso_FCTransportProblem* TransportProblemAdapter::getPaso_FCTransportProblem() const
57 {
58 return m_transport_problem.get();
59 }
60
61
62 void TransportProblemAdapter::setToSolution(escript::Data& out, escript::Data& source, const double dt, const boost::python::dict& options) const
63 {
64 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
65 Paso_Options paso_options;
66 SystemMatrixAdapter::dictToPasoOptions(&paso_options,options);
67 if ( out.getDataPointSize() != getBlockSize()) {
68 throw FinleyAdapterException("solve : block size of solution does not match block size of transport problems.");
69 } else if ( source.getDataPointSize() != getBlockSize()) {
70 throw FinleyAdapterException("solve : block size of source term does not match block size of transport problems.");
71 } else if ( out.getFunctionSpace() != getFunctionSpace()) {
72 throw FinleyAdapterException("solve : function spaces of solution and of transport problem don't match.");
73 } else if (source.getFunctionSpace() != getFunctionSpace()) {
74 throw FinleyAdapterException("solve : function spaces of source term and of transport problem don't match.");
75 } else if (dt<=0.) {
76 throw FinleyAdapterException("solve : time increment dt needs to be positive.");
77 }
78 out.expand();
79 source.expand();
80 double* out_dp=out.getSampleData(0);
81 double* source_dp=source.getSampleData(0);
82 Paso_SolverFCT_solve(transp,out_dp,dt,source_dp,&paso_options);
83 checkPasoError();
84 }
85
86 void TransportProblemAdapter::resetTransport() const
87 {
88 Paso_FCTransportProblem* transp = getPaso_FCTransportProblem();
89 Paso_FCTransportProblem_reset(transp);
90 checkPasoError();
91 }
92 void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r) const
93 {
94 if ( q.getDataPointSize() != getBlockSize()) {
95 throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint mask.");
96 } else if ( q.getFunctionSpace() != getFunctionSpace()) {
97 throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
98 } else if ( r.getDataPointSize() != getBlockSize()) {
99 throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint values.");
100 } else if ( r.getFunctionSpace() != getFunctionSpace()) {
101 throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint values don't match.");
102 } else if ( source.getDataPointSize() != getBlockSize()) {
103 throw FinleyAdapterException("copyConstraint : block size does not match the number of components of source.");
104 } else if ( source.getFunctionSpace() != getFunctionSpace()) {
105 throw FinleyAdapterException("copyConstraint : function spaces of transport problem and source don't match.");
106 }
107 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
108
109 /* r2=r where q>0, 0 elsewhere */
110 escript::Data r2(0.,q.getDataPointShape(),q.getFunctionSpace());
111 r2.copyWithMask(r,q);
112
113 /* source-=transp->mass_matrix*r2 */
114 double* r2_dp=r2.getSampleData(0);
115 double* source_dp=source.getSampleData(0);
116 r2.expand();
117 source.expand();
118 Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
119 checkPasoError();
120
121 /* insert 0 rows and cols into transport matrix */
122 q.expand();
123 double* q_dp=q.getSampleData(0);
124 Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
125 checkPasoError();
126
127 /* insert 0 rows amd 1 in main diagonal into mass matrix */
128 Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
129 checkPasoError();
130
131 source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
132 }
133
134 void TransportProblemAdapter::copyInitialValue(escript::Data& u) const
135 {
136 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
137 if ( u.getDataPointSize() != getBlockSize()) {
138 throw FinleyAdapterException("copyInitialValue : block size of solution does not match block size of transport problems.");
139 } else if ( u.getFunctionSpace() != getFunctionSpace()) {
140 throw FinleyAdapterException("copyInitialValue : function spaces of solution and of transport problem don't match.");
141 }
142 u.expand();
143 double* u_dp=u.getSampleData(0);
144 Paso_FCTransportProblem_checkinSolution( transp,u_dp);
145 checkPasoError();
146 }
147
148 double TransportProblemAdapter::getSafeTimeStepSize() const
149 {
150 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
151 double dt=Paso_FCTransportProblem_getSafeTimeStepSize(transp);
152 checkPasoError();
153 return dt;
154 }
155
156 double TransportProblemAdapter::getUnlimitedTimeStepSize() const
157 {
158 return LARGE_POSITIVE_FLOAT;
159 }
160
161
162
163 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26