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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3823 - (hide annotations)
Mon Feb 13 08:09:32 2012 UTC (7 years, 10 months ago) by gross
File MIME type: text/plain
File size: 5524 byte(s)
some fixes for openmp/mpi compilation and linearNC can now handel constraints.
1 gross 2987
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 caltinay 3642 /* Paso: TransportProblem (see Paso_TransportSolver_solve) */
18 gross 2987
19     /**************************************************************/
20    
21 caltinay 3642 /* Author: l.gross@uq.edu.au */
22 gross 2987
23     /**************************************************************/
24    
25    
26     #include "Transport.h"
27     #include "PasoUtil.h"
28    
29    
30     /**************************************************************/
31    
32 caltinay 3642 /* free all used memory */
33 gross 2987
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 jfenwick 3259 Esys_MPIInfo_free(in->mpi_info);
42 gross 3822 MEMFREE(in->constraint_mask);
43 gross 3005 MEMFREE(in->reactive_matrix);
44     MEMFREE(in->main_diagonal_mass_matrix);
45 gross 2987 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 gross 3793 Paso_TransportProblem* Paso_TransportProblem_alloc(Paso_SystemMatrixPattern *pattern, const int block_size)
60 gross 2987 {
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 gross 3793
65 gross 2987 out=MEMALLOC(1,Paso_TransportProblem);
66 jfenwick 3259 if (Esys_checkPtr(out)) return NULL;
67 gross 3793 out->reference_counter=0;
68     out->dt_max_R=LARGE_POSITIVE_FLOAT;
69     out->dt_max_T=LARGE_POSITIVE_FLOAT;
70 gross 2987 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 gross 3822 out->constraint_mask=NULL;
75 jfenwick 3259 out->mpi_info=Esys_MPIInfo_getReference(pattern->mpi_info);
76 gross 3005
77 gross 2987 out->lumped_mass_matrix=NULL;
78     out->main_diagonal_low_order_transport_matrix=NULL;
79 gross 3005 out->reactive_matrix=NULL;
80     out->main_diagonal_mass_matrix=NULL;
81 gross 2987
82 jfenwick 3259 if (Esys_noError()) {
83 gross 2987 n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix);
84 gross 3822 out->constraint_mask=MEMALLOC(n,double); /* ? */
85 gross 3793 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 gross 2987
90 gross 3793
91 gross 3822 if ( ! (Esys_checkPtr(out->constraint_mask) ||
92 jfenwick 3259 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 gross 2987 {
95 gross 3005 #pragma omp parallel for schedule(static) private(i)
96 gross 2987 for (i = 0; i < n; ++i) {
97     out->lumped_mass_matrix[i]=0.;
98     out->main_diagonal_low_order_transport_matrix[i]=0.;
99 gross 3822 out->constraint_mask[i]=0.;
100 gross 2987 }
101 gross 3005 }
102 gross 2987 }
103 jfenwick 3259 if (Esys_noError()) {
104 gross 2987 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 gross 3822 const dim_t n = Paso_SystemMatrix_getTotalNumRows(in->transport_matrix);
115 gross 2987 Paso_SystemMatrix_setValues(in->transport_matrix, 0.);
116     Paso_SystemMatrix_setValues(in->mass_matrix, 0.);
117     Paso_solve_free(in->iteration_matrix);
118 gross 3822 Paso_zeroes( n, in->constraint_mask);
119 gross 2987 in->valid_matrices=FALSE;
120     }
121    
122    
123 jfenwick 3259 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 gross 2987 {
125     return MATRIX_FORMAT_DEFAULT + MATRIX_FORMAT_BLK1;
126     }
127    
128 gross 3822 void Paso_TransportProblem_setUpConstraint(Paso_TransportProblem* fctp, const double* q)
129 gross 2987 {
130 gross 3005 dim_t i;
131 gross 3822 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
132 gross 2987 if ( fctp->valid_matrices ) {
133 caltinay 3642 Esys_setError(VALUE_ERROR, "Paso_TransportProblem_setUpConstraint: Cannot insert a constraint into a valid system.");
134 gross 2987 return;
135     }
136    
137 gross 3823 #pragma omp parallel for schedule(static) private(i)
138 gross 2987 for (i=0;i<n;++i) {
139     if (q[i]>0) {
140 gross 3822 fctp->constraint_mask[i]=1;
141 gross 2987 } else {
142 gross 3822 fctp->constraint_mask[i]=0;
143 gross 2987 }
144     }
145 gross 3822 fctp->valid_matrices=FALSE;
146 gross 2987 }
147    
148     void Paso_TransportProblem_insertConstraint(Paso_TransportProblem* fctp, const double* r, double* source)
149     {
150 gross 3822 dim_t i;
151     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
152 gross 2987
153 gross 3026 #pragma omp parallel for schedule(static) private(i)
154 gross 3822 for (i=0;i<n;++i) {
155     if (fctp->constraint_mask[i] > 0) source[i] = r[i];
156     }
157 gross 2987 }
158 gross 3822

  ViewVC Help
Powered by ViewVC 1.1.26