/[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 1513 - (show annotations)
Tue Apr 15 08:47:57 2008 UTC (11 years, 2 months ago) by gross
File size: 6181 byte(s)
pragma ivdep removed. icc produced wrong code.
1 /*******************************************************
2 *
3 * Copyright 2007 by University of Queensland
4 *
5 * http://esscc.uq.edu.au
6 * Primary Business: Queensland, Australia
7 * Licensed under the Open Software License version 3.0
8 * http://www.opensource.org/licenses/osl-3.0.php
9 *
10 *******************************************************/
11
12 #ifdef PASO_MPI
13 #include <mpi.h>
14 #endif
15 #include "TransportProblemAdapter.h"
16 #include "SystemMatrixAdapter.h"
17
18 using namespace std;
19
20 namespace finley {
21
22 struct null_deleter
23 {
24 void operator()(void const *ptr) const
25 {
26 }
27 };
28
29
30 TransportProblemAdapter::TransportProblemAdapter()
31 {
32 throw FinleyAdapterException("Error - Illegal to generate default TransportProblemAdapter.");
33 }
34
35
36 TransportProblemAdapter::TransportProblemAdapter(Paso_FCTransportProblem* transport_problem,
37 const double theta,
38 const int block_size,
39 const escript::FunctionSpace& functionspace):
40 AbstractTransportProblem(theta, block_size, functionspace)
41 {
42 m_transport_problem.reset(transport_problem,null_deleter());
43 }
44
45 TransportProblemAdapter::~TransportProblemAdapter()
46 {
47 if (m_transport_problem.unique()) {
48 Paso_FCTransportProblem* transp=m_transport_problem.get();
49 Paso_FCTransportProblem_free(transp);
50 }
51 }
52
53 Paso_FCTransportProblem* TransportProblemAdapter::getPaso_FCTransportProblem() const
54 {
55 return m_transport_problem.get();
56 }
57
58
59 void TransportProblemAdapter::setToSolution(escript::Data& out, escript::Data& source, const double dt, const boost::python::dict& options) const
60 {
61 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
62 Paso_Options paso_options;
63 SystemMatrixAdapter::dictToPasoOptions(&paso_options,options);
64 if ( out.getDataPointSize() != getBlockSize()) {
65 throw FinleyAdapterException("solve : block size of solution does not match block size of transport problems.");
66 } else if ( source.getDataPointSize() != getBlockSize()) {
67 throw FinleyAdapterException("solve : block size of source term does not match block size of transport problems.");
68 } else if ( out.getFunctionSpace() != getFunctionSpace()) {
69 throw FinleyAdapterException("solve : function spaces of solution and of transport problem don't match.");
70 } else if (source.getFunctionSpace() != getFunctionSpace()) {
71 throw FinleyAdapterException("solve : function spaces of source term and of transport problem don't match.");
72 } else if (dt<=0.) {
73 throw FinleyAdapterException("solve : time increment dt needs to be positive.");
74 }
75 out.expand();
76 source.expand();
77 double* out_dp=out.getSampleData(0);
78 double* source_dp=source.getSampleData(0);
79 Paso_SolverFCT_solve(transp,out_dp,dt,source_dp,&paso_options);
80 checkPasoError();
81 }
82
83 void TransportProblemAdapter::resetTransport() const
84 {
85 Paso_FCTransportProblem* transp = getPaso_FCTransportProblem();
86 Paso_FCTransportProblem_reset(transp);
87 checkPasoError();
88 }
89 void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r) const
90 {
91 if ( q.getDataPointSize() != getBlockSize()) {
92 throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint mask.");
93 } else if ( q.getFunctionSpace() != getFunctionSpace()) {
94 throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
95 } else if ( r.getDataPointSize() != getBlockSize()) {
96 throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint values.");
97 } else if ( r.getFunctionSpace() != getFunctionSpace()) {
98 throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint values don't match.");
99 } else if ( source.getDataPointSize() != getBlockSize()) {
100 throw FinleyAdapterException("copyConstraint : block size does not match the number of components of source.");
101 } else if ( source.getFunctionSpace() != getFunctionSpace()) {
102 throw FinleyAdapterException("copyConstraint : function spaces of transport problem and source don't match.");
103 }
104 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
105
106 /* r2=r where q>0, 0 elsewhere */
107 escript::Data r2(0.,q.getDataPointShape(),q.getFunctionSpace());
108 r2.copyWithMask(r,q);
109
110 /* source-=transp->mass_matrix*r2 */
111 double* r2_dp=r2.getSampleData(0);
112 double* source_dp=source.getSampleData(0);
113 r2.expand();
114 source.expand();
115 Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
116 checkPasoError();
117
118 /* insert 0 rows and cols into transport matrix */
119 q.expand();
120 double* q_dp=q.getSampleData(0);
121 Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
122 checkPasoError();
123
124 /* insert 0 rows amd 1 in main diagonal into mass matrix */
125 Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
126 checkPasoError();
127
128 source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
129 }
130
131 void TransportProblemAdapter::copyInitialValue(escript::Data& u) const
132 {
133 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
134 if ( u.getDataPointSize() != getBlockSize()) {
135 throw FinleyAdapterException("copyInitialValue : block size of solution does not match block size of transport problems.");
136 } else if ( u.getFunctionSpace() != getFunctionSpace()) {
137 throw FinleyAdapterException("copyInitialValue : function spaces of solution and of transport problem don't match.");
138 }
139 u.expand();
140 double* u_dp=u.getSampleData(0);
141 Paso_FCTransportProblem_checkinSolution( transp,u_dp);
142 checkPasoError();
143 }
144
145 double TransportProblemAdapter::getSafeTimeStepSize() const
146 {
147 Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
148 double dt=Paso_FCTransportProblem_getSafeTimeStepSize(transp);
149 checkPasoError();
150 return dt;
151 }
152
153
154
155 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26