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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2315 - (show annotations)
Wed Mar 18 00:38:48 2009 UTC (10 years, 7 months ago) by gross
File MIME type: text/plain
File size: 8676 byte(s)
fixes for MPI
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 /**************************************************************/
16
17 /* Paso: FCTransportProblem */
18
19 /**************************************************************/
20
21 /* Author: l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25
26 #include "Paso.h"
27 #include "SolverFCT.h"
28 #include "PasoUtil.h"
29
30
31 /**************************************************************/
32
33 /* free all memory used by */
34
35 void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
36 if (in!=NULL) {
37 in->reference_counter--;
38 if (in->reference_counter<=0) {
39 Paso_SystemMatrix_free(in->transport_matrix);
40 Paso_SystemMatrix_free(in->mass_matrix);
41 Paso_SystemMatrix_free(in->iteration_matrix);
42 Paso_MPIInfo_free(in->mpi_info);
43 Paso_Coupler_free(in->u_coupler);
44 MEMFREE(in->u);
45 MEMFREE(in->constraint_weights);
46 MEMFREE(in->main_iptr);
47 MEMFREE(in->lumped_mass_matrix);
48 MEMFREE(in->main_diagonal_low_order_transport_matrix);
49 MEMFREE(in);
50 }
51 }
52 }
53
54 Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
55 if (in!=NULL) {
56 ++(in->reference_counter);
57 }
58 return in;
59 }
60
61 Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
62 return in->transport_matrix;
63 }
64 Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
65 return in->mass_matrix;
66 }
67
68 double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
69 return in->lumped_mass_matrix;
70 }
71
72 dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
73 return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
74 }
75
76 Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size)
77 {
78 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
79 Paso_FCTransportProblem* out=NULL;
80 Paso_SystemMatrixPattern *transport_pattern;
81 dim_t n,i;
82 index_t iptr,iptr_main;
83
84 if ((theta<0.) || (theta >1.)) {
85 Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
86 return NULL;
87 }
88
89 out=MEMALLOC(1,Paso_FCTransportProblem);
90 if (Paso_checkPtr(out)) return NULL;
91 out->reference_counter=0;
92 out->theta=theta;
93 out->dt_max=LARGE_POSITIVE_FLOAT;
94 out->constraint_factor=sqrt(LARGE_POSITIVE_FLOAT);
95 if (out->theta < 1.) {
96 out->dt_factor=MIN(1./(1.-out->theta),DT_FACTOR_MAX);
97 } else {
98 out->dt_factor=DT_FACTOR_MAX;
99 }
100 out->valid_matrices=FALSE;
101 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
102 out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
103 out->iteration_matrix=NULL;
104 out->u=NULL;
105 out->constraint_weights=NULL;
106 out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
107 out->main_iptr=NULL;
108 out->lumped_mass_matrix=NULL;
109 out->main_diagonal_low_order_transport_matrix=NULL;
110
111 if (Paso_noError()) {
112 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
113 transport_pattern=out->transport_matrix->pattern;
114 out->u=MEMALLOC(n,double);
115 out->constraint_weights=MEMALLOC(n,double);
116 out->main_iptr=MEMALLOC(n,index_t);
117 out->lumped_mass_matrix=MEMALLOC(n,double);
118 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
119 out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);
120
121 if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) || Paso_checkPtr(out->constraint_weights) ||
122 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError() ) {
123
124 #pragma omp parallel
125 {
126 #pragma omp for schedule(static) private(i)
127 for (i = 0; i < n; ++i) {
128 out->lumped_mass_matrix[i]=0.;
129 out->main_diagonal_low_order_transport_matrix[i]=0.;
130 out->u[i]=0.;
131 }
132 /* identify the main diagonals */
133 #pragma omp for schedule(static) private(i,iptr,iptr_main)
134 for (i = 0; i < n; ++i) {
135 iptr_main=transport_pattern->mainPattern->ptr[0]-1;
136 for (iptr=transport_pattern->mainPattern->ptr[i];iptr<transport_pattern->mainPattern->ptr[i+1]; iptr++) {
137 if (transport_pattern->mainPattern->index[iptr]==i) {
138 iptr_main=iptr;
139 break;
140 }
141 }
142 out->main_iptr[i]=iptr_main;
143 if (iptr_main==transport_pattern->mainPattern->ptr[0]-1)
144 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
145 }
146
147 }
148 }
149
150 }
151 if (Paso_noError()) {
152 out->reference_counter=1;
153 return out;
154 } else {
155 Paso_FCTransportProblem_free(out);
156 return NULL;
157 }
158 }
159
160 void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
161 dim_t i, n;
162 double local_u_min,u_min;
163 n=Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
164 u_min=LARGE_POSITIVE_FLOAT;
165 #pragma omp parallel private(local_u_min)
166 {
167 local_u_min=0.;
168 #pragma omp for schedule(static) private(i)
169 for (i=0;i<n;++i) {
170 in->u[i]=u[i];
171 local_u_min=MIN(local_u_min,u[i]);
172 }
173 #pragma omp critical
174 {
175 u_min=MIN(u_min,local_u_min);
176 }
177 }
178 #ifdef PASO_MPI
179 local_u_min=u_min;
180 MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
181 #endif
182 if (u_min<0) {
183 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_checkinSolution: initial guess must be non-negative.");
184 }
185 }
186 dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
187 {
188 return in->transport_matrix->row_block_size;
189 }
190
191 Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
192 {
193 return in->transport_matrix->pattern->col_connector;
194 }
195
196 index_t Paso_FCTransportProblem_getTypeId(const index_t solver,const index_t preconditioner, const index_t package,const bool_t symmetry, Paso_MPIInfo *mpi_info)
197 {
198 return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
199 }
200
201 void Paso_FCTransportProblem_setUpConstraint(Paso_FCTransportProblem* fctp, const double* q, const double factor)
202 {
203 dim_t i, n;
204 register double m, rtmp;
205 double factor2= fctp->dt_factor * factor;
206 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
207
208 if ( fctp->valid_matrices ) {
209 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: you must not insert a constraint is a valid system.");
210 return;
211 }
212 if (factor<=0) {
213 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_insertConstraint: constraint_factor needs to be positive.");
214 return;
215 }
216
217
218 #pragma omp for schedule(static) private(i,m,rtmp)
219 for (i=0;i<n;++i) {
220 if (q[i]>0) {
221 m=fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]];
222 rtmp=factor2 * (m == 0 ? 1 : m);
223 fctp->constraint_weights[i]=rtmp;
224 fctp->mass_matrix->mainBlock->val[fctp->main_iptr[i]]=m+rtmp;
225 } else {
226 fctp->constraint_weights[i]=0;
227 }
228 }
229 fctp->constraint_factor=factor;
230 }
231
232 void Paso_FCTransportProblem_insertConstraint(Paso_FCTransportProblem* fctp, const double* r, double* source)
233 {
234 dim_t i, n;
235 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
236
237 #pragma omp for schedule(static) private(i)
238 for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
239 }

  ViewVC Help
Powered by ViewVC 1.1.26