src/CPPAdapter/TransportProblemAdapter.cpp


/*******************************************************
*
* Copyright (c) 2003-2009 by University of Queensland
* Earth Systems Science Computational Center (ESSCC)
* http://www.uq.edu.au/esscc
*
* Primary Business: Queensland, Australia
* Licensed under the Open Software License version 3.0
* http://www.opensource.org/licenses/osl-3.0.php
*
*******************************************************/


#ifdef PASO_MPI
#include <mpi.h>
#endif
#include "TransportProblemAdapter.h" 
#include "SystemMatrixAdapter.h" 

using namespace std;

namespace finley {

struct null_deleter
{
  void operator()(void const *ptr) const
  {
  }
};


TransportProblemAdapter::TransportProblemAdapter()
{
   throw FinleyAdapterException("Error - Illegal to generate default TransportProblemAdapter.");
}


TransportProblemAdapter::TransportProblemAdapter(Paso_FCTransportProblem* transport_problem,
                                                 const double theta,
                                                 const int block_size,
                                                 const escript::FunctionSpace& functionspace):
AbstractTransportProblem(theta, block_size, functionspace)
{
    m_transport_problem.reset(transport_problem,null_deleter());
}

TransportProblemAdapter::~TransportProblemAdapter()
{ 
    if (m_transport_problem.unique()) {
        Paso_FCTransportProblem* transp=m_transport_problem.get();
        Paso_FCTransportProblem_free(transp);
    }
}

Paso_FCTransportProblem* TransportProblemAdapter::getPaso_FCTransportProblem() const 
{
   return m_transport_problem.get();
}


void TransportProblemAdapter::setToSolution(escript::Data& out, escript::Data& source, const double dt, boost::python::object& options) const
{
    Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
    Paso_Options paso_options;
    SystemMatrixAdapter::escriptToPasoOptions(&paso_options,options);
    options.attr("resetDiagnostics")();
    if ( out.getDataPointSize()  != getBlockSize()) {
     throw FinleyAdapterException("solve : block size of solution does not match block size of transport problems.");
    } else if ( source.getDataPointSize() != getBlockSize()) {
     throw FinleyAdapterException("solve : block size of source term does not match block size of transport problems.");
    } else if ( out.getFunctionSpace()  != getFunctionSpace()) {
     throw FinleyAdapterException("solve : function spaces of solution and of transport problem don't match.");
    } else if (source.getFunctionSpace() != getFunctionSpace()) {
     throw FinleyAdapterException("solve : function spaces of source term and of transport problem don't match.");
    } else if (dt<=0.) {
     throw FinleyAdapterException("solve : time increment dt needs to be positive.");
    }
    out.expand();
    source.expand();
    out.requireWrite();
    source.requireWrite();
    double* out_dp=out.getSampleDataRW(0);
    double* source_dp=source.getSampleDataRW(0);
    Paso_SolverFCT_solve(transp,out_dp,dt,source_dp,&paso_options);
    SystemMatrixAdapter::pasoToEscriptOptions(&paso_options,options);
    checkPasoError();
}

void TransportProblemAdapter::resetTransport() const
{
   Paso_FCTransportProblem* transp = getPaso_FCTransportProblem();
   Paso_FCTransportProblem_reset(transp);
   checkPasoError();
}
void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r, const double factor) const
{
    if ( q.getDataPointSize()  != getBlockSize()) {
     throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint mask.");
    } else if ( q.getFunctionSpace()  != getFunctionSpace()) {
     throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
    } else if ( r.getDataPointSize()  != getBlockSize()) {
     throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint values.");
    } else if ( r.getFunctionSpace()  != getFunctionSpace()) {
     throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint values don't match.");
    } else if ( source.getDataPointSize()  != getBlockSize()) {
     throw FinleyAdapterException("copyConstraint : block size does not match the number of components of source.");
    } else if ( source.getFunctionSpace()  != getFunctionSpace()) {
     throw FinleyAdapterException("copyConstraint : function spaces of transport problem and source don't match.");
    }
    Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();

    /* r2=r where q>0, 0 elsewhere */
    escript::Data r2(0.,q.getDataPointShape(),q.getFunctionSpace());
    r2.copyWithMask(r,q);

    /* source-=transp->mass_matrix*r2 */
    r2.expand();
    source.expand();
    q.expand();
    r2.requireWrite();
    source.requireWrite();
    q.requireWrite();
    double* r2_dp=r2.getSampleDataRW(0);
    double* source_dp=source.getSampleDataRW(0);
    double* q_dp=q.getSampleDataRW(0);

    if (false) {
       cout << "v1\n";
       Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
       checkPasoError();

       /* insert 0 rows into transport matrix */
       Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
       checkPasoError();

       /* insert 0 rows amd 1 in main diagonal into mass matrix */
       Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
       checkPasoError();

       source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
   } else {
       Paso_FCTransportProblem_setUpConstraint(transp, q_dp, factor);
       checkPasoError();
       Paso_FCTransportProblem_insertConstraint(transp,r2_dp, source_dp);
       checkPasoError();
   }
}

void TransportProblemAdapter::copyInitialValue(escript::Data& u) const
{
    Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
    if ( u.getDataPointSize()  != getBlockSize()) {
     throw FinleyAdapterException("copyInitialValue : block size of solution does not match block size of transport problems.");
    } else if ( u.getFunctionSpace()  != getFunctionSpace()) {
     throw FinleyAdapterException("copyInitialValue : function spaces of solution and of transport problem don't match.");
    }
    u.expand();
    u.requireWrite();
    double* u_dp=u.getSampleDataRW(0);
    Paso_FCTransportProblem_checkinSolution( transp,u_dp);
    checkPasoError();
}

double TransportProblemAdapter::getSafeTimeStepSize() const
{
    Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
    double dt=Paso_FCTransportProblem_getSafeTimeStepSize(transp);
    checkPasoError();
    return dt;
}

double TransportProblemAdapter::getUnlimitedTimeStepSize() const
{
    return LARGE_POSITIVE_FLOAT;
}


}  // end of namespace
1
2	/*******************************************************
3	*
4	* Copyright (c) 2003-2009 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, boost::python::object& options) const
63	{
64	Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
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* source_dp=source.getSampleDataRW(0);
85	Paso_SolverFCT_solve(transp,out_dp,dt,source_dp,&paso_options);
86	SystemMatrixAdapter::pasoToEscriptOptions(&paso_options,options);
87	checkPasoError();
88	}
89
90	void TransportProblemAdapter::resetTransport() const
91	{
92	Paso_FCTransportProblem* transp = getPaso_FCTransportProblem();
93	Paso_FCTransportProblem_reset(transp);
94	checkPasoError();
95	}
96	void TransportProblemAdapter::copyConstraint(escript::Data& source, escript::Data& q, escript::Data& r, const double factor) const
97	{
98	if ( q.getDataPointSize() != getBlockSize()) {
99	throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint mask.");
100	} else if ( q.getFunctionSpace() != getFunctionSpace()) {
101	throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint mask don't match.");
102	} else if ( r.getDataPointSize() != getBlockSize()) {
103	throw FinleyAdapterException("copyConstraint : block size does not match the number of components of constraint values.");
104	} else if ( r.getFunctionSpace() != getFunctionSpace()) {
105	throw FinleyAdapterException("copyConstraint : function spaces of transport problem and constraint values don't match.");
106	} else if ( source.getDataPointSize() != getBlockSize()) {
107	throw FinleyAdapterException("copyConstraint : block size does not match the number of components of source.");
108	} else if ( source.getFunctionSpace() != getFunctionSpace()) {
109	throw FinleyAdapterException("copyConstraint : function spaces of transport problem and source don't match.");
110	}
111	Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
112
113	/* r2=r where q>0, 0 elsewhere */
114	escript::Data r2(0.,q.getDataPointShape(),q.getFunctionSpace());
115	r2.copyWithMask(r,q);
116
117	/* source-=transp->mass_matrixr2 /
118	r2.expand();
119	source.expand();
120	q.expand();
121	r2.requireWrite();
122	source.requireWrite();
123	q.requireWrite();
124	double* r2_dp=r2.getSampleDataRW(0);
125	double* source_dp=source.getSampleDataRW(0);
126	double* q_dp=q.getSampleDataRW(0);
127
128	if (false) {
129	cout << "v1\n";
130	Paso_SystemMatrix_MatrixVector(-1., transp->mass_matrix, r2_dp, 1., source_dp);
131	checkPasoError();
132
133	/* insert 0 rows into transport matrix */
134	Paso_SystemMatrix_nullifyRows(transp->transport_matrix,q_dp, 0.);
135	checkPasoError();
136
137	/* insert 0 rows amd 1 in main diagonal into mass matrix */
138	Paso_SystemMatrix_nullifyRowsAndCols(transp->mass_matrix,q_dp,q_dp,1.);
139	checkPasoError();
140
141	source.copyWithMask(escript::Data(0.,q.getDataPointShape(),q.getFunctionSpace()),q);
142	} else {
143	Paso_FCTransportProblem_setUpConstraint(transp, q_dp, factor);
144	checkPasoError();
145	Paso_FCTransportProblem_insertConstraint(transp,r2_dp, source_dp);
146	checkPasoError();
147	}
148	}
149
150	void TransportProblemAdapter::copyInitialValue(escript::Data& u) const
151	{
152	Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
153	if ( u.getDataPointSize() != getBlockSize()) {
154	throw FinleyAdapterException("copyInitialValue : block size of solution does not match block size of transport problems.");
155	} else if ( u.getFunctionSpace() != getFunctionSpace()) {
156	throw FinleyAdapterException("copyInitialValue : function spaces of solution and of transport problem don't match.");
157	}
158	u.expand();
159	u.requireWrite();
160	double* u_dp=u.getSampleDataRW(0);
161	Paso_FCTransportProblem_checkinSolution( transp,u_dp);
162	checkPasoError();
163	}
164
165	double TransportProblemAdapter::getSafeTimeStepSize() const
166	{
167	Paso_FCTransportProblem* transp=getPaso_FCTransportProblem();
168	double dt=Paso_FCTransportProblem_getSafeTimeStepSize(transp);
169	checkPasoError();
170	return dt;
171	}
172
173	double TransportProblemAdapter::getUnlimitedTimeStepSize() const
174	{
175	return LARGE_POSITIVE_FLOAT;
176	}
177
178
179
180	} // end of namespace
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision