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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1401 - (show annotations)
Fri Jan 25 04:31:18 2008 UTC (12 years, 8 months ago) by gross
File MIME type: text/plain
File size: 5647 byte(s)
rewrite antidiffusion calculation to avoid coloring for OPENMP parallelization
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
38 Paso_SystemMatrix_free(in->transport_matrix);
39 Paso_SystemMatrix_free(in->flux_matrix);
40 Paso_MPIInfo_free(in->mpi_info);
41
42 MEMFREE(in->u);
43 MEMFREE(in->lumped_mass_matrix);
44 MEMFREE(in->row_sum_flux_matrix);
45 MEMFREE(in->transport_matrix_diagonal);
46 MEMFREE(in->r_p);
47 MEMFREE(in->r_n);
48 MEMFREE(in->main_iptr);
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
65 Paso_SystemMatrix* Paso_FCTransportProblem_borrowFluxMatrix(Paso_FCTransportProblem* in) {
66 return in->flux_matrix;
67 }
68
69 double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
70 return in->lumped_mass_matrix;
71 }
72
73 dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
74 return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
75 }
76
77 Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, double dt_max, Paso_SystemMatrixPattern *pattern, int block_size
78
79
80 ) {
81 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
82 Paso_FCTransportProblem* out=NULL;
83 dim_t n,i;
84 index_t iptr,iptr_main,k;
85
86 if ((theta<0.) || (theta >1.)) {
87 Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
88 return NULL;
89 }
90
91 out=MEMALLOC(1,Paso_FCTransportProblem);
92 if (Paso_checkPtr(out)) return NULL;
93
94 out->theta=theta;
95 out->dt_max=dt_max;
96 out->valid_matrices=FALSE;
97 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
98 Paso_SystemMatrix_allocBuffer(out->transport_matrix);
99 out->flux_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
100 out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
101
102 out->main_iptr=NULL;
103 out->lumped_mass_matrix=NULL;
104 out->row_sum_flux_matrix=NULL;
105 out->transport_matrix_diagonal=NULL;
106 out->r_p=NULL;
107 out->r_n=NULL;
108
109 if (Paso_noError()) {
110 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
111
112 out->r_p=MEMALLOC(n,double);
113 out->r_n=MEMALLOC(n,double);
114 out->main_iptr=MEMALLOC(n,index_t);
115 out->lumped_mass_matrix=MEMALLOC(n,double);
116 out->row_sum_flux_matrix=MEMALLOC(n,double);
117 out->transport_matrix_diagonal=MEMALLOC(n,double);
118 out->u=MEMALLOC(n,double);
119
120 if ( ! (Paso_checkPtr(out->r_p) || Paso_checkPtr(out->r_n) || Paso_checkPtr(out->main_iptr) ||
121 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->transport_matrix_diagonal) || Paso_checkPtr(out->row_sum_flux_matrix) || Paso_checkPtr(out->u)) ) {
122
123 #pragma omp parallel for schedule(static) private(i)
124 for (i = 0; i < n; ++i) {
125 out->lumped_mass_matrix[i]=0.;
126 out->row_sum_flux_matrix[i]=0.;
127 out->u[i]=0.;
128 out->r_p[i]=0.;
129 out->r_n[i]=0.;
130 }
131 /* identify the main diagonals */
132 #pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k)
133 for (i = 0; i < n; ++i) {
134 iptr_main=pattern->mainPattern->ptr[0]-1;
135 for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
136 if (pattern->mainPattern->index[iptr]==i) {
137 iptr_main=iptr;
138 break;
139 }
140 }
141 out->main_iptr[i]=iptr_main;
142 if (iptr_main==pattern->mainPattern->ptr[0]-1)
143 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
144 }
145
146 }
147
148 }
149 if (Paso_noError()) {
150 return out;
151 } else {
152 Paso_FCTransportProblem_free(out);
153 return NULL;
154 }
155 }
156
157 void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
158 dim_t i, n;
159
160 n=Paso_FCTransportProblem_getTotalNumRows(in);
161 #pragma omp parallel for schedule(static) private(i)
162 for (i = 0; i < n; ++i) {
163 in->u[i]=u[i];
164 }
165 }

  ViewVC Help
Powered by ViewVC 1.1.26