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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1552 - (show annotations)
Thu May 8 08:52:41 2008 UTC (11 years, 9 months ago) by gross
File MIME type: text/plain
File size: 5847 byte(s)
some changes to make the implementatiopn of a upwind MPI version easier
1 /* $Id:$ */
2
3 /*******************************************************
4 *
5 * Copyright 2007 by University of Queensland
6 *
7 * http://esscc.uq.edu.au
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 /* Paso: FCTransportProblem */
17
18 /**************************************************************/
19
20 /* Author: l.gross@uq.edu.au */
21
22 /**************************************************************/
23
24
25 #include "Paso.h"
26 #include "SolverFCT.h"
27 #include "PasoUtil.h"
28
29
30 /**************************************************************/
31
32 /* free all memory used by */
33 void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
34 if (in!=NULL) {
35 in->reference_counter--;
36 if (in->reference_counter<=0) {
37 Paso_SystemMatrix_free(in->transport_matrix);
38 Paso_SystemMatrix_free(in->mass_matrix);
39 Paso_SystemMatrix_free(in->iteration_matrix);
40 Paso_MPIInfo_free(in->mpi_info);
41 MEMFREE(in->u);
42 MEMFREE(in->main_iptr);
43 MEMFREE(in->lumped_mass_matrix);
44 MEMFREE(in->main_diagonal_low_order_transport_matrix);
45 MEMFREE(in);
46 }
47 }
48 }
49
50 Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
51 if (in!=NULL) {
52 ++(in->reference_counter);
53 }
54 return in;
55 }
56
57 Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
58 return in->transport_matrix;
59 }
60 Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
61 return in->mass_matrix;
62 }
63
64 double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
65 return in->lumped_mass_matrix;
66 }
67
68 dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
69 return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
70 }
71
72 Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size
73
74
75 ) {
76 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
77 Paso_FCTransportProblem* out=NULL;
78 dim_t n,i;
79 index_t iptr,iptr_main,k;
80
81 if ((theta<0.) || (theta >1.)) {
82 Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
83 return NULL;
84 }
85
86 out=MEMALLOC(1,Paso_FCTransportProblem);
87 if (Paso_checkPtr(out)) return NULL;
88
89
90 out->theta=theta;
91 out->dt_max=LARGE_POSITIVE_FLOAT;
92 out->valid_matrices=FALSE;
93 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
94 out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
95 out->iteration_matrix=NULL;
96
97 out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
98 out->u=NULL;
99 out->u_min=0.;
100 out->main_iptr=NULL;
101 out->lumped_mass_matrix=NULL;
102 out->main_diagonal_low_order_transport_matrix=NULL;
103
104 if (Paso_noError()) {
105 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
106
107 out->u=MEMALLOC(n,double);
108 out->main_iptr=MEMALLOC(n,index_t);
109 out->lumped_mass_matrix=MEMALLOC(n,double);
110 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
111
112 if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
113 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) ) {
114
115 #pragma omp parallel for schedule(static) private(i)
116 for (i = 0; i < n; ++i) {
117 out->lumped_mass_matrix[i]=0.;
118 out->main_diagonal_low_order_transport_matrix[i]=0.;
119 out->u[i]=0.;
120 }
121 /* identify the main diagonals */
122 #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
123 for (i = 0; i < n; ++i) {
124 iptr_main=pattern->mainPattern->ptr[0]-1;
125 for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
126 if (pattern->mainPattern->index[iptr]==i) {
127 iptr_main=iptr;
128 break;
129 }
130 }
131 out->main_iptr[i]=iptr_main;
132 if (iptr_main==pattern->mainPattern->ptr[0]-1)
133 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
134 }
135
136 }
137
138 }
139 if (Paso_noError()) {
140 return out;
141 } else {
142 Paso_FCTransportProblem_free(out);
143 return NULL;
144 }
145 }
146
147 void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
148 dim_t i, n;
149 double local_u_min,u_min;
150
151 n=Paso_FCTransportProblem_getTotalNumRows(in);
152 u_min=LARGE_POSITIVE_FLOAT;
153 #pragma omp parallel private(local_u_min)
154 {
155 local_u_min=0.;
156 #pragma omp for schedule(static) private(i)
157 for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]);
158 #pragma omp critical
159 {
160 u_min=MIN(u_min,local_u_min);
161 }
162 }
163 #ifdef PASO_MPI
164 local_u_min=u_min;
165 MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
166 #endif
167 in->u_min=u_min;
168 #pragma omp parallel for schedule(static) private(i)
169 for (i=0;i<n;++i) {
170 in->u[i]=u[i]-u_min;
171 }
172 }

  ViewVC Help
Powered by ViewVC 1.1.26