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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2987 - (show annotations)
Tue Mar 16 01:32:43 2010 UTC (9 years, 10 months ago) by gross
File MIME type: text/plain
File size: 8162 byte(s)
FCT solver rewritten
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 Paso_MPIInfo_free(in->mpi_info);
42 Paso_Coupler_free(in->u_coupler);
43 MEMFREE(in->constraint_weights);
44 MEMFREE(in->main_iptr);
45 MEMFREE(in->lumped_mass_matrix);
46 MEMFREE(in->main_diagonal_low_order_transport_matrix);
47 MEMFREE(in);
48 }
49 }
50 }
51
52 Paso_TransportProblem* Paso_TransportProblem_getReference(Paso_TransportProblem* in) {
53 if (in!=NULL) {
54 ++(in->reference_counter);
55 }
56 return in;
57 }
58
59 Paso_SystemMatrix* Paso_TransportProblem_borrowTransportMatrix(Paso_TransportProblem* in) {
60 return in->transport_matrix;
61 }
62 Paso_SystemMatrix* Paso_TransportProblem_borrowMassMatrix(Paso_TransportProblem* in) {
63 return in->mass_matrix;
64 }
65
66 double* Paso_TransportProblem_borrowLumpedMassMatrix(Paso_TransportProblem* in) {
67 return in->lumped_mass_matrix;
68 }
69
70 dim_t Paso_TransportProblem_getTotalNumRows(Paso_TransportProblem* in) {
71 return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
72 }
73
74 Paso_TransportProblem* Paso_TransportProblem_alloc(bool_t useBackwardEuler, Paso_SystemMatrixPattern *pattern, int block_size)
75 {
76 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
77 Paso_TransportProblem* out=NULL;
78 Paso_SystemMatrixPattern *transport_pattern;
79 dim_t n,i;
80 int fail=0;
81 index_t iptr,iptr_main;
82
83 out=MEMALLOC(1,Paso_TransportProblem);
84 if (Paso_checkPtr(out)) return NULL;
85 out->reference_counter=0;
86 out->useBackwardEuler=useBackwardEuler;
87 out->dt_max=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=Paso_MPIInfo_getReference(pattern->mpi_info);
102 out->main_iptr=NULL;
103 out->lumped_mass_matrix=NULL;
104 out->main_diagonal_low_order_transport_matrix=NULL;
105
106 if (Paso_noError()) {
107 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
108 transport_pattern=out->transport_matrix->pattern;
109 out->constraint_weights=MEMALLOC(n,double);
110 out->main_iptr=MEMALLOC(n,index_t);
111 out->lumped_mass_matrix=MEMALLOC(n,double);
112 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
113 out->u_coupler=Paso_Coupler_alloc(Paso_TransportProblem_borrowConnector(out),block_size);
114
115 if ( ! (Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->constraint_weights) ||
116 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError() )
117 {
118 fail=0;
119 #pragma omp parallel
120 {
121 #pragma omp for schedule(static) private(i)
122 for (i = 0; i < n; ++i) {
123 out->lumped_mass_matrix[i]=0.;
124 out->main_diagonal_low_order_transport_matrix[i]=0.;
125 }
126 /* identify the main diagonals */
127 #pragma omp for schedule(static) private(i,iptr,iptr_main)
128 for (i = 0; i < n; ++i) {
129 iptr_main=transport_pattern->mainPattern->ptr[0]-1;
130 for (iptr=transport_pattern->mainPattern->ptr[i];iptr<transport_pattern->mainPattern->ptr[i+1]; iptr++) {
131 if (transport_pattern->mainPattern->index[iptr]==i) {
132 iptr_main=iptr;
133 break;
134 }
135 }
136 out->main_iptr[i]=iptr_main;
137 if (iptr_main==transport_pattern->mainPattern->ptr[0]-1) fail=1;
138 }
139
140 }
141 #ifdef PASO_MPI
142 {
143 int fail_loc = fail;
144 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, out->mpi_info->comm);
145 }
146 #endif
147 if (fail>0) Paso_setError(VALUE_ERROR, "Paso_TransportProblem_alloc: no main diagonal");
148 }
149 }
150 if (Paso_noError()) {
151 out->reference_counter=1;
152 return out;
153 } else {
154 Paso_TransportProblem_free(out);
155 return NULL;
156 }
157 }
158
159 void Paso_TransportProblem_reset(Paso_TransportProblem* in)
160 {
161 Paso_SystemMatrix_setValues(in->transport_matrix, 0.);
162 Paso_SystemMatrix_setValues(in->mass_matrix, 0.);
163 Paso_solve_free(in->iteration_matrix);
164 in->valid_matrices=FALSE;
165 }
166
167 dim_t Paso_TransportProblem_getBlockSize(const Paso_TransportProblem* in)
168 {
169 return in->transport_matrix->row_block_size;
170 }
171
172 Paso_Connector* Paso_TransportProblem_borrowConnector(const Paso_TransportProblem* in)
173 {
174 return in->transport_matrix->pattern->col_connector;
175 }
176
177 index_t Paso_TransportProblem_getTypeId(const index_t solver,const index_t preconditioner, const index_t package,const bool_t symmetry, Paso_MPIInfo *mpi_info)
178 {
179 return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
180 }
181
182 /****************** REVISE ****************************/
183 void Paso_TransportProblem_setUpConstraint(Paso_TransportProblem* fctp, const double* q, const double factor)
184 {
185 dim_t i, n;
186 register double m, rtmp;
187 double factor2= fctp->dt_factor * factor;
188 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
189
190 if ( fctp->valid_matrices ) {
191 Paso_setError(VALUE_ERROR, "Paso_TransportProblem_insertConstraint: you must not insert a constraint is a valid system.");
192 return;
193 }
194 if (factor<=0) {
195 Paso_setError(VALUE_ERROR, "Paso_TransportProblem_insertConstraint: constraint_factor needs to be positive.");
196 return;
197 }
198
199
200 #pragma omp for schedule(static) private(i,m,rtmp)
201 for (i=0;i<n;++i) {
202 if (q[i]>0) {
203 m=fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]];
204 rtmp=factor2 * (m == 0 ? 1 : m);
205 fctp->constraint_weights[i]=rtmp;
206 fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]]=m+rtmp;
207 } else {
208 fctp->constraint_weights[i]=0;
209 }
210 }
211 fctp->constraint_factor=factor;
212 }
213
214 void Paso_TransportProblem_insertConstraint(Paso_TransportProblem* fctp, const double* r, double* source)
215 {
216 dim_t i, n;
217 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
218
219 #pragma omp for schedule(static) private(i)
220 for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
221 }
222 /*****************************************************/

  ViewVC Help
Powered by ViewVC 1.1.26