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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 5 months ago) by ksteube
File MIME type: text/plain
File size: 6539 byte(s)
Copyright updated in all files

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->main_iptr);
46 MEMFREE(in->lumped_mass_matrix);
47 MEMFREE(in->main_diagonal_low_order_transport_matrix);
48 MEMFREE(in);
49 }
50 }
51 }
52
53 Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
54 if (in!=NULL) {
55 ++(in->reference_counter);
56 }
57 return in;
58 }
59
60 Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
61 return in->transport_matrix;
62 }
63 Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
64 return in->mass_matrix;
65 }
66
67 double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
68 return in->lumped_mass_matrix;
69 }
70
71 dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
72 return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
73 }
74
75 Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size)
76 {
77 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
78 Paso_FCTransportProblem* out=NULL;
79 dim_t n,i;
80 index_t iptr,iptr_main;
81
82 if ((theta<0.) || (theta >1.)) {
83 Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
84 return NULL;
85 }
86
87 out=MEMALLOC(1,Paso_FCTransportProblem);
88 if (Paso_checkPtr(out)) return NULL;
89 out->reference_counter=0;
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 out->u=NULL;
97 out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
98 out->main_iptr=NULL;
99 out->lumped_mass_matrix=NULL;
100 out->main_diagonal_low_order_transport_matrix=NULL;
101
102 if (Paso_noError()) {
103 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
104
105 out->u=MEMALLOC(n,double);
106 out->main_iptr=MEMALLOC(n,index_t);
107 out->lumped_mass_matrix=MEMALLOC(n,double);
108 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
109 out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);
110
111 if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
112 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError() ) {
113
114 #pragma omp parallel
115 {
116 #pragma omp for schedule(static) private(i)
117 for (i = 0; i < n; ++i) {
118 out->lumped_mass_matrix[i]=0.;
119 out->main_diagonal_low_order_transport_matrix[i]=0.;
120 out->u[i]=0.;
121 }
122 /* identify the main diagonals */
123 #pragma omp for schedule(static) private(i,iptr,iptr_main)
124 for (i = 0; i < n; ++i) {
125 iptr_main=pattern->mainPattern->ptr[0]-1;
126 for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) {
127 if (pattern->mainPattern->index[iptr]==i) {
128 iptr_main=iptr;
129 break;
130 }
131 }
132 out->main_iptr[i]=iptr_main;
133 if (iptr_main==pattern->mainPattern->ptr[0]-1)
134 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal");
135 }
136
137 }
138 }
139
140 }
141 if (Paso_noError()) {
142 out->reference_counter=1;
143 return out;
144 } else {
145 Paso_FCTransportProblem_free(out);
146 return NULL;
147 }
148 }
149
150 void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
151 dim_t i, n;
152 double local_u_min,u_min;
153 n=Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
154 u_min=LARGE_POSITIVE_FLOAT;
155 #pragma omp parallel private(local_u_min)
156 {
157 local_u_min=0.;
158 #pragma omp for schedule(static) private(i)
159 for (i=0;i<n;++i) {
160 in->u[i]=u[i];
161 local_u_min=MIN(local_u_min,u[i]);
162 }
163 #pragma omp critical
164 {
165 u_min=MIN(u_min,local_u_min);
166 }
167 }
168 #ifdef PASO_MPI
169 local_u_min=u_min;
170 MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
171 #endif
172 if (u_min<0) {
173 Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_checkinSolution: initial guess must be non-negative.");
174 }
175 }
176 dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
177 {
178 return in->transport_matrix->row_block_size;
179 }
180
181 Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
182 {
183 return in->transport_matrix->pattern->col_connector;
184 }

  ViewVC Help
Powered by ViewVC 1.1.26