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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3793 - (show annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 6 months ago) by gross
File MIME type: text/plain
File size: 6233 byte(s)
new implementation of FCT solver with some modifications to the python interface
1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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: TransportProblem (see Paso_TransportSolver_solve) */
18
19 /**************************************************************/
20
21 /* Author: l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25
26 #include "Transport.h"
27 #include "PasoUtil.h"
28
29
30 /**************************************************************/
31
32 /* free all used memory */
33
34 void Paso_TransportProblem_free(Paso_TransportProblem* in) {
35 if (in!=NULL) {
36 in->reference_counter--;
37 if (in->reference_counter<=0) {
38 Paso_SystemMatrix_free(in->transport_matrix);
39 Paso_SystemMatrix_free(in->mass_matrix);
40 Paso_SystemMatrix_free(in->iteration_matrix);
41 Esys_MPIInfo_free(in->mpi_info);
42 MEMFREE(in->constraint_weights);
43 MEMFREE(in->reactive_matrix);
44 MEMFREE(in->main_diagonal_mass_matrix);
45 MEMFREE(in->lumped_mass_matrix);
46 MEMFREE(in->main_diagonal_low_order_transport_matrix);
47 MEMFREE(in);
48 }
49 }
50 }
51
52 Paso_TransportProblem* Paso_TransportProblem_getReference(Paso_TransportProblem* in) {
53 if (in!=NULL) {
54 ++(in->reference_counter);
55 }
56 return in;
57 }
58
59 Paso_TransportProblem* Paso_TransportProblem_alloc(Paso_SystemMatrixPattern *pattern, const int block_size)
60 {
61 Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */
62 Paso_TransportProblem* out=NULL;
63 dim_t n,i;
64
65 out=MEMALLOC(1,Paso_TransportProblem);
66 if (Esys_checkPtr(out)) return NULL;
67 out->reference_counter=0;
68 out->dt_max_R=LARGE_POSITIVE_FLOAT;
69 out->dt_max_T=LARGE_POSITIVE_FLOAT;
70
71
72
73 /****************** REVISE ****************************/
74 out->constraint_factor=sqrt(LARGE_POSITIVE_FLOAT);
75
76 out->dt_factor=2.;
77 /*****************************************************/
78 out->valid_matrices=FALSE;
79 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size,FALSE);
80 out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size,FALSE);
81 out->iteration_matrix=NULL;
82 out->constraint_weights=NULL;
83 out->mpi_info=Esys_MPIInfo_getReference(pattern->mpi_info);
84
85 out->lumped_mass_matrix=NULL;
86 out->main_diagonal_low_order_transport_matrix=NULL;
87 out->reactive_matrix=NULL;
88 out->main_diagonal_mass_matrix=NULL;
89
90 if (Esys_noError()) {
91 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
92 out->constraint_weights=MEMALLOC(n,double); /* ? */
93 out->lumped_mass_matrix=MEMALLOC(n,double); /* ? */
94 out->reactive_matrix=MEMALLOC(n,double); /* ? */
95 out->main_diagonal_mass_matrix=MEMALLOC(n,double); /* ? */
96 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double); /* ? */
97
98
99 if ( ! (Esys_checkPtr(out->constraint_weights) ||
100 Esys_checkPtr(out->reactive_matrix) || Esys_checkPtr(out->main_diagonal_mass_matrix) ||
101 Esys_checkPtr(out->lumped_mass_matrix) || Esys_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Esys_noError() )
102 {
103 #pragma omp parallel for schedule(static) private(i)
104 for (i = 0; i < n; ++i) {
105 out->lumped_mass_matrix[i]=0.;
106 out->main_diagonal_low_order_transport_matrix[i]=0.;
107 }
108 }
109 }
110 if (Esys_noError()) {
111 out->reference_counter=1;
112 return out;
113 } else {
114 Paso_TransportProblem_free(out);
115 return NULL;
116 }
117 }
118
119 void Paso_TransportProblem_reset(Paso_TransportProblem* in)
120 {
121 Paso_SystemMatrix_setValues(in->transport_matrix, 0.);
122 Paso_SystemMatrix_setValues(in->mass_matrix, 0.);
123 Paso_solve_free(in->iteration_matrix);
124 in->valid_matrices=FALSE;
125 }
126
127
128 index_t Paso_TransportProblem_getTypeId(const index_t solver,const index_t preconditioner, const index_t package,const bool_t symmetry, Esys_MPIInfo *mpi_info)
129 {
130 return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
131 }
132
133 /****************** REVISE ****************************/
134 void Paso_TransportProblem_setUpConstraint(Paso_TransportProblem* fctp, const double* q, const double factor)
135 {
136 dim_t i;
137 register double m, rtmp;
138 double factor2= fctp->dt_factor * factor;
139 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
140 const index_t* main_iptr=Paso_SparseMatrix_borrowMainDiagonalPointer(fctp->mass_matrix->mainBlock);
141
142 if ( fctp->valid_matrices ) {
143 Esys_setError(VALUE_ERROR, "Paso_TransportProblem_setUpConstraint: Cannot insert a constraint into a valid system.");
144 return;
145 }
146 if (factor<=0) {
147 Esys_setError(VALUE_ERROR, "Paso_TransportProblem_setUpConstraint: constraint_factor needs to be positive.");
148 return;
149 }
150
151
152 #pragma omp parallel for schedule(static) private(i,m,rtmp)
153 for (i=0;i<n;++i) {
154 if (q[i]>0) {
155 m=fctp->mass_matrix->mainBlock->val[main_iptr[i]];
156 rtmp=factor2 * (m == 0 ? 1 : m);
157 fctp->constraint_weights[i]=rtmp;
158 fctp->mass_matrix->mainBlock->val[main_iptr[i]]=m+rtmp;
159 } else {
160 fctp->constraint_weights[i]=0;
161 }
162 }
163 fctp->constraint_factor=factor;
164 }
165
166 void Paso_TransportProblem_insertConstraint(Paso_TransportProblem* fctp, const double* r, double* source)
167 {
168 dim_t i, n;
169 n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
170
171 #pragma omp parallel for schedule(static) private(i)
172 for (i=0;i<n;++i) source[i]+=fctp->constraint_weights[i] * r[i];
173 }
174 /*****************************************************/

  ViewVC Help
Powered by ViewVC 1.1.26