/[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 2197 - (show annotations)
Thu Jan 8 05:49:16 2009 UTC (10 years, 3 months ago) by gross
File size: 6601 byte(s)
modifications to the way cosntraints are handeled in the Transport problem


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 double factor) 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 r2.expand();
115 source.expand();
116 q.expand();
117 double* r2_dp=r2.getSampleData(0);
118 double* source_dp=source.getSampleData(0);
119 double* q_dp=q.getSampleData(0);
120
121 if (false) {
122 cout << "v1\n";
123 Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
124 checkPasoError();
125
126 /* insert 0 rows into transport matrix */
127 Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
128 checkPasoError();
129
130 /* insert 0 rows amd 1 in main diagonal into mass matrix */
131 Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
132 checkPasoError();
133
134 source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
135 } else {
136 Paso_FCTransportProblem_setUpConstraint(transp, q_dp, factor);
137 checkPasoError();
138 Paso_FCTransportProblem_insertConstraint(transp,r2_dp, source_dp);
139 checkPasoError();
140 }
141 }
142
143 void TransportProblemAdapter::copyInitialValue(escript::Data& u) const
144 {
145 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
146 if ( u.getDataPointSize() != getBlockSize()) {
147 throw FinleyAdapterException("copyInitialValue : block size of solution does not match block size of transport problems.");
148 } else if ( u.getFunctionSpace() != getFunctionSpace()) {
149 throw FinleyAdapterException("copyInitialValue : function spaces of solution and of transport problem don't match.");
150 }
151 u.expand();
152 double* u_dp=u.getSampleData(0);
153 Paso_FCTransportProblem_checkinSolution( transp,u_dp);
154 checkPasoError();
155 }
156
157 double TransportProblemAdapter::getSafeTimeStepSize() const
158 {
159 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
160 double dt=Paso_FCTransportProblem_getSafeTimeStepSize(transp);
161 checkPasoError();
162 return dt;
163 }
164
165 double TransportProblemAdapter::getUnlimitedTimeStepSize() const
166 {
167 return LARGE_POSITIVE_FLOAT;
168 }
169
170
171
172 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26