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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1401 - (hide annotations)
Fri Jan 25 04:31:18 2008 UTC (12 years, 1 month ago) by gross
File MIME type: text/plain
File size: 5647 byte(s)
rewrite antidiffusion calculation to avoid coloring for OPENMP parallelization
1 gross 1361 /* $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 gross 1362 /* Paso: FCTransportProblem */
17 gross 1361
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 gross 1362 /**************************************************************/
31 gross 1361
32 gross 1362 /* 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 gross 1361
38 gross 1362 Paso_SystemMatrix_free(in->transport_matrix);
39     Paso_SystemMatrix_free(in->flux_matrix);
40     Paso_MPIInfo_free(in->mpi_info);
41 gross 1361
42 gross 1364 MEMFREE(in->u);
43 gross 1362 MEMFREE(in->lumped_mass_matrix);
44 gross 1363 MEMFREE(in->row_sum_flux_matrix);
45 gross 1370 MEMFREE(in->transport_matrix_diagonal);
46 gross 1401 MEMFREE(in->r_p);
47     MEMFREE(in->r_n);
48 gross 1362 MEMFREE(in->main_iptr);
49     MEMFREE(in);
50     }
51     }
52     }
53 gross 1361
54 gross 1362 Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
55 gross 1361 if (in!=NULL) {
56 gross 1362 ++(in->reference_counter);
57 gross 1361 }
58 gross 1367 return in;
59 gross 1362 }
60    
61     Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
62     return in->transport_matrix;
63 gross 1361 }
64    
65 gross 1362 Paso_SystemMatrix* Paso_FCTransportProblem_borrowFluxMatrix(Paso_FCTransportProblem* in) {
66     return in->flux_matrix;
67     }
68 gross 1361
69 gross 1362 double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
70     return in->lumped_mass_matrix;
71     }
72 gross 1361
73 gross 1362 dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
74     return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
75     }
76 gross 1361
77 gross 1366 Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, double dt_max, Paso_SystemMatrixPattern *pattern, int block_size
78 gross 1361
79 gross 1362
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 gross 1361 return NULL;
89 gross 1362 }
90 gross 1361
91 gross 1362 out=MEMALLOC(1,Paso_FCTransportProblem);
92     if (Paso_checkPtr(out)) return NULL;
93 gross 1361
94 gross 1362 out->theta=theta;
95 gross 1370 out->dt_max=dt_max;
96 gross 1363 out->valid_matrices=FALSE;
97 gross 1362 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 gross 1361
102 gross 1362 out->main_iptr=NULL;
103     out->lumped_mass_matrix=NULL;
104 gross 1363 out->row_sum_flux_matrix=NULL;
105 gross 1370 out->transport_matrix_diagonal=NULL;
106 gross 1401 out->r_p=NULL;
107     out->r_n=NULL;
108 gross 1361
109 gross 1362 if (Paso_noError()) {
110     n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
111    
112 gross 1401 out->r_p=MEMALLOC(n,double);
113     out->r_n=MEMALLOC(n,double);
114 gross 1362 out->main_iptr=MEMALLOC(n,index_t);
115     out->lumped_mass_matrix=MEMALLOC(n,double);
116 gross 1363 out->row_sum_flux_matrix=MEMALLOC(n,double);
117 gross 1370 out->transport_matrix_diagonal=MEMALLOC(n,double);
118 gross 1364 out->u=MEMALLOC(n,double);
119 gross 1362
120 gross 1401 if ( ! (Paso_checkPtr(out->r_p) || Paso_checkPtr(out->r_n) || Paso_checkPtr(out->main_iptr) ||
121 gross 1370 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->transport_matrix_diagonal) || Paso_checkPtr(out->row_sum_flux_matrix) || Paso_checkPtr(out->u)) ) {
122 gross 1362
123 gross 1367 #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 gross 1401 out->r_p[i]=0.;
129     out->r_n[i]=0.;
130 gross 1367 }
131 gross 1362 /* 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 gross 1361 }
149     if (Paso_noError()) {
150     return out;
151     } else {
152 gross 1362 Paso_FCTransportProblem_free(out);
153 gross 1361 return NULL;
154     }
155     }
156 gross 1364
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