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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1804 - (hide annotations)
Wed Sep 24 07:52:19 2008 UTC (11 years, 5 months ago) by gross
File MIME type: text/plain
File size: 6526 byte(s)
a robister version of the upwinding scheme
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 gross 1804
34 gross 1362 void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) {
35     if (in!=NULL) {
36     in->reference_counter--;
37     if (in->reference_counter<=0) {
38     Paso_SystemMatrix_free(in->transport_matrix);
39 gross 1407 Paso_SystemMatrix_free(in->mass_matrix);
40     Paso_SystemMatrix_free(in->iteration_matrix);
41 gross 1362 Paso_MPIInfo_free(in->mpi_info);
42 gross 1804 Paso_Coupler_free(in->u_coupler);
43 gross 1364 MEMFREE(in->u);
44 gross 1407 MEMFREE(in->main_iptr);
45 gross 1362 MEMFREE(in->lumped_mass_matrix);
46 gross 1407 MEMFREE(in->main_diagonal_low_order_transport_matrix);
47 gross 1362 MEMFREE(in);
48     }
49     }
50     }
51 gross 1361
52 gross 1362 Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) {
53 gross 1361 if (in!=NULL) {
54 gross 1362 ++(in->reference_counter);
55 gross 1361 }
56 gross 1367 return in;
57 gross 1362 }
58    
59     Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) {
60     return in->transport_matrix;
61 gross 1361 }
62 gross 1407 Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) {
63     return in->mass_matrix;
64 gross 1362 }
65 gross 1361
66 gross 1362 double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) {
67     return in->lumped_mass_matrix;
68     }
69 gross 1361
70 gross 1362 dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) {
71     return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
72     }
73 gross 1361
74 gross 1804 Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size)
75     {
76 gross 1362 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 phornby 1628 index_t iptr,iptr_main;
80 gross 1362
81     if ((theta<0.) || (theta >1.)) {
82     Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1.");
83 gross 1361 return NULL;
84 gross 1362 }
85 gross 1361
86 gross 1362 out=MEMALLOC(1,Paso_FCTransportProblem);
87     if (Paso_checkPtr(out)) return NULL;
88 gross 1804 out->reference_counter=0;
89 gross 1362 out->theta=theta;
90 gross 1407 out->dt_max=LARGE_POSITIVE_FLOAT;
91 gross 1363 out->valid_matrices=FALSE;
92 gross 1362 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
93 gross 1407 out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size);
94     out->iteration_matrix=NULL;
95     out->u=NULL;
96 gross 1562 out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
97 gross 1362 out->main_iptr=NULL;
98     out->lumped_mass_matrix=NULL;
99 gross 1407 out->main_diagonal_low_order_transport_matrix=NULL;
100 gross 1361
101 gross 1362 if (Paso_noError()) {
102     n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
103    
104 gross 1407 out->u=MEMALLOC(n,double);
105 gross 1362 out->main_iptr=MEMALLOC(n,index_t);
106     out->lumped_mass_matrix=MEMALLOC(n,double);
107 gross 1407 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double);
108 gross 1804 out->u_coupler=Paso_Coupler_alloc(Paso_FCTransportProblem_borrowConnector(out),block_size);
109 gross 1362
110 gross 1407 if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) ||
111 gross 1804 Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Paso_noError() ) {
112 gross 1362
113 gross 1562 #pragma omp parallel
114     {
115     #pragma omp 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 phornby 1628 #pragma omp for schedule(static) private(i,iptr,iptr_main)
123 gross 1562 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 gross 1367 }
137 gross 1362 }
138    
139 gross 1361 }
140     if (Paso_noError()) {
141 gross 1804 out->reference_counter=1;
142 gross 1361 return out;
143     } else {
144 gross 1362 Paso_FCTransportProblem_free(out);
145 gross 1361 return NULL;
146     }
147     }
148 gross 1364
149     void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) {
150     dim_t i, n;
151 gross 1407 double local_u_min,u_min;
152 gross 1804 n=Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
153 gross 1407 u_min=LARGE_POSITIVE_FLOAT;
154     #pragma omp parallel private(local_u_min)
155     {
156     local_u_min=0.;
157     #pragma omp for schedule(static) private(i)
158 gross 1804 for (i=0;i<n;++i) {
159     in->u[i]=u[i];
160     local_u_min=MIN(local_u_min,u[i]);
161     }
162 gross 1407 #pragma omp critical
163     {
164     u_min=MIN(u_min,local_u_min);
165     }
166 gross 1364 }
167 gross 1407 #ifdef PASO_MPI
168     local_u_min=u_min;
169 ksteube 1415 MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm);
170 gross 1407 #endif
171 gross 1804 if (u_min<0) {
172     Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_checkinSolution: initial guess must be non-negative.");
173 gross 1410 }
174 gross 1364 }
175 gross 1562 dim_t Paso_FCTransportProblem_getBlockSize(const Paso_FCTransportProblem* in)
176     {
177     return in->transport_matrix->row_block_size;
178     }
179    
180     Paso_Connector* Paso_FCTransportProblem_borrowConnector(const Paso_FCTransportProblem* in)
181     {
182     return in->transport_matrix->pattern->col_connector;
183     }

  ViewVC Help
Powered by ViewVC 1.1.26