/[escript]/branches/symbolic_from_3470/paso/src/Transport.c
ViewVC logotype

Contents of /branches/symbolic_from_3470/paso/src/Transport.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3868 - (show annotations)
Thu Mar 15 06:07:08 2012 UTC (7 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 5524 byte(s)
Update to latest trunk

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_mask);
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 out->valid_matrices=FALSE;
71 out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size,FALSE);
72 out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size,FALSE);
73 out->iteration_matrix=NULL;
74 out->constraint_mask=NULL;
75 out->mpi_info=Esys_MPIInfo_getReference(pattern->mpi_info);
76
77 out->lumped_mass_matrix=NULL;
78 out->main_diagonal_low_order_transport_matrix=NULL;
79 out->reactive_matrix=NULL;
80 out->main_diagonal_mass_matrix=NULL;
81
82 if (Esys_noError()) {
83 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
84 out->constraint_mask=MEMALLOC(n,double); /* ? */
85 out->lumped_mass_matrix=MEMALLOC(n,double); /* ? */
86 out->reactive_matrix=MEMALLOC(n,double); /* ? */
87 out->main_diagonal_mass_matrix=MEMALLOC(n,double); /* ? */
88 out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double); /* ? */
89
90
91 if ( ! (Esys_checkPtr(out->constraint_mask) ||
92 Esys_checkPtr(out->reactive_matrix) || Esys_checkPtr(out->main_diagonal_mass_matrix) ||
93 Esys_checkPtr(out->lumped_mass_matrix) || Esys_checkPtr(out->main_diagonal_low_order_transport_matrix)) && Esys_noError() )
94 {
95 #pragma omp parallel for schedule(static) private(i)
96 for (i = 0; i < n; ++i) {
97 out->lumped_mass_matrix[i]=0.;
98 out->main_diagonal_low_order_transport_matrix[i]=0.;
99 out->constraint_mask[i]=0.;
100 }
101 }
102 }
103 if (Esys_noError()) {
104 out->reference_counter=1;
105 return out;
106 } else {
107 Paso_TransportProblem_free(out);
108 return NULL;
109 }
110 }
111
112 void Paso_TransportProblem_reset(Paso_TransportProblem* in)
113 {
114 const dim_t n = Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
115 Paso_SystemMatrix_setValues(in->transport_matrix, 0.);
116 Paso_SystemMatrix_setValues(in->mass_matrix, 0.);
117 Paso_solve_free(in->iteration_matrix);
118 Paso_zeroes( n, in->constraint_mask);
119 in->valid_matrices=FALSE;
120 }
121
122
123 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)
124 {
125 return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
126 }
127
128 void Paso_TransportProblem_setUpConstraint(Paso_TransportProblem* fctp, const double* q)
129 {
130 dim_t i;
131 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
132 if ( fctp->valid_matrices ) {
133 Esys_setError(VALUE_ERROR, "Paso_TransportProblem_setUpConstraint: Cannot insert a constraint into a valid system.");
134 return;
135 }
136
137 #pragma omp parallel for schedule(static) private(i)
138 for (i=0;i<n;++i) {
139 if (q[i]>0) {
140 fctp->constraint_mask[i]=1;
141 } else {
142 fctp->constraint_mask[i]=0;
143 }
144 }
145 fctp->valid_matrices=FALSE;
146 }
147
148 void Paso_TransportProblem_insertConstraint(Paso_TransportProblem* fctp, const double* r, double* source)
149 {
150 dim_t i;
151 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
152
153 #pragma omp parallel for schedule(static) private(i)
154 for (i=0;i<n;++i) {
155 if (fctp->constraint_mask[i] > 0) source[i] = r[i];
156 }
157 }
158

  ViewVC Help
Powered by ViewVC 1.1.26