/[escript]/trunk/paso/src/Transport.c
ViewVC logotype

Contents of /trunk/paso/src/Transport.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 7367 byte(s)
Merging dudley and scons updates from branches

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 /**************************************************************/
16
17 /* Paso: TransportProblem (see Paso_TransportSolver_solve) */
18
19 /**************************************************************/
20
21 /* Author: l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25
26 #include "Transport.h"
27 #include "PasoUtil.h"
28
29
30 /**************************************************************/
31
32 /* free all memory used by */
33
34 void Paso_TransportProblem_free(Paso_TransportProblem* in) {
35 if (in!=NULL) {
36 in->reference_counter--;
37 if (in->reference_counter<=0) {
38 Paso_SystemMatrix_free(in->transport_matrix);
39 Paso_SystemMatrix_free(in->mass_matrix);
40 Paso_SystemMatrix_free(in->iteration_matrix);
41 Esys_MPIInfo_free(in->mpi_info);
42 Paso_Coupler_free(in->u_coupler);
43 MEMFREE(in->constraint_weights);
44 MEMFREE(in->reactive_matrix);
45 MEMFREE(in->main_diagonal_mass_matrix);
46 MEMFREE(in->lumped_mass_matrix);
47 MEMFREE(in->main_diagonal_low_order_transport_matrix);
48 MEMFREE(in);
49 }
50 }
51 }
52
53 Paso_TransportProblem* Paso_TransportProblem_getReference(Paso_TransportProblem* in) {
54 if (in!=NULL) {
55 ++(in->reference_counter);
56 }
57 return in;
58 }
59
60 Paso_SystemMatrix* Paso_TransportProblem_borrowTransportMatrix(Paso_TransportProblem* in) {
61 return in->transport_matrix;
62 }
63 Paso_SystemMatrix* Paso_TransportProblem_borrowMassMatrix(Paso_TransportProblem* in) {
64 return in->mass_matrix;
65 }
66
67 double* Paso_TransportProblem_borrowLumpedMassMatrix(Paso_TransportProblem* in) {
68 return in->lumped_mass_matrix;
69 }
70
71 dim_t Paso_TransportProblem_getTotalNumRows(Paso_TransportProblem* in) {
72 return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
73 }
74
75 Paso_TransportProblem* Paso_TransportProblem_alloc(bool_t useBackwardEuler, Paso_SystemMatrixPattern *pattern, int block_size)
76 {
77 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
78 Paso_TransportProblem* out=NULL;
79 Paso_SystemMatrixPattern *transport_pattern;
80 dim_t n,i;
81
82 out=MEMALLOC(1,Paso_TransportProblem);
83 if (Esys_checkPtr(out)) return NULL;
84 out->reference_counter=0;
85 out->useBackwardEuler=useBackwardEuler;
86 out->dt_max=LARGE_POSITIVE_FLOAT;
87 out->dt_failed=LARGE_POSITIVE_FLOAT;
88 /****************** REVISE ****************************/
89 out->constraint_factor=sqrt(LARGE_POSITIVE_FLOAT);
90 if (out->useBackwardEuler) {
91 out->dt_factor=DT_FACTOR_MAX;
92 } else {
93 out->dt_factor=2.;
94 }
95 /*****************************************************/
96 out->valid_matrices=FALSE;
97 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size,FALSE);
98 out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size,FALSE);
99 out->iteration_matrix=NULL;
100 out->constraint_weights=NULL;
101 out->mpi_info=Esys_MPIInfo_getReference(pattern->mpi_info);
102
103 out->lumped_mass_matrix=NULL;
104 out->main_diagonal_low_order_transport_matrix=NULL;
105 out->reactive_matrix=NULL;
106 out->main_diagonal_mass_matrix=NULL;
107
108 if (Esys_noError()) {
109 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
110 transport_pattern=out->transport_matrix->pattern;
111 out->constraint_weights=MEMALLOC(n,double);
112 out->lumped_mass_matrix=MEMALLOC(n,double);
113 out->reactive_matrix=MEMALLOC(n,double);;
114 out->main_diagonal_mass_matrix=MEMALLOC(n,double);
115 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
116 out->u_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(out),block_size);
117
118 if ( ! (Esys_checkPtr(out->constraint_weights) ||
119 Esys_checkPtr(out->reactive_matrix) || Esys_checkPtr(out->main_diagonal_mass_matrix) ||
120 Esys_checkPtr(out->lumped_mass_matrix) || Esys_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Esys_noError() )
121 {
122 #pragma omp parallel for schedule(static) private(i)
123 for (i = 0; i < n; ++i) {
124 out->lumped_mass_matrix[i]=0.;
125 out->main_diagonal_low_order_transport_matrix[i]=0.;
126 }
127 }
128 }
129 if (Esys_noError()) {
130 out->reference_counter=1;
131 return out;
132 } else {
133 Paso_TransportProblem_free(out);
134 return NULL;
135 }
136 }
137
138 void Paso_TransportProblem_reset(Paso_TransportProblem* in)
139 {
140 Paso_SystemMatrix_setValues(in->transport_matrix, 0.);
141 Paso_SystemMatrix_setValues(in->mass_matrix, 0.);
142 Paso_solve_free(in->iteration_matrix);
143 in->valid_matrices=FALSE;
144 }
145
146 dim_t Paso_TransportProblem_getBlockSize(const Paso_TransportProblem* in)
147 {
148 return in->transport_matrix->row_block_size;
149 }
150
151 Paso_Connector* Paso_TransportProblem_borrowConnector(const Paso_TransportProblem* in)
152 {
153 return in->transport_matrix->pattern->col_connector;
154 }
155
156 index_t Paso_TransportProblem_getTypeId(const index_t solver,const index_t preconditioner, const index_t package,const bool_t symmetry, Esys_MPIInfo *mpi_info)
157 {
158 return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
159 }
160
161 /****************** REVISE ****************************/
162 void Paso_TransportProblem_setUpConstraint(Paso_TransportProblem* fctp, const double* q, const double factor)
163 {
164 dim_t i;
165 register double m, rtmp;
166 double factor2= fctp->dt_factor * factor;
167 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
168 const index_t* main_iptr=Paso_SparseMatrix_borrowMainDiagonalPointer(fctp->mass_matrix->mainBlock);
169
170 if ( fctp->valid_matrices ) {
171 Esys_setError(VALUE_ERROR, "Paso_TransportProblem_insertConstraint: you must not insert a constraint is a valid system.");
172 return;
173 }
174 if (factor<=0) {
175 Esys_setError(VALUE_ERROR, "Paso_TransportProblem_insertConstraint: constraint_factor needs to be positive.");
176 return;
177 }
178
179
180 #pragma omp parallel for schedule(static) private(i,m,rtmp)
181 for (i=0;i<n;++i) {
182 if (q[i]>0) {
183 m=fctp->mass_matrix->mainBlock->val[main_iptr[i]];
184 rtmp=factor2 * (m == 0 ? 1 : m);
185 fctp->constraint_weights[i]=rtmp;
186 fctp->mass_matrix->mainBlock->val[main_iptr[i]]=m+rtmp;
187 } else {
188 fctp->constraint_weights[i]=0;
189 }
190 }
191 fctp->constraint_factor=factor;
192 }
193
194 void Paso_TransportProblem_insertConstraint(Paso_TransportProblem* fctp, const double* r, double* source)
195 {
196 dim_t i, n;
197 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
198
199 #pragma omp parallel for schedule(static) private(i)
200 for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
201 }
202 /*****************************************************/

  ViewVC Help
Powered by ViewVC 1.1.26