/[escript]/trunk/paso/src/FluxLimiter.cpp
ViewVC logotype

Contents of /trunk/paso/src/FluxLimiter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4934 - (show annotations)
Tue May 13 00:28:11 2014 UTC (5 years, 4 months ago) by jfenwick
File size: 8987 byte(s)
This commit is brought to you by the number 4934 and the tool "meld".
Merge of partially complete split world code from branch.



1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************/
19
20 /* Paso: CFT: fluxlimiter for a transport problem
21 *
22 ****************************************************************************/
23
24 /* Author: l.gross@uq.edu.au */
25
26 /****************************************************************************/
27
28 #include "FluxLimiter.h"
29
30 namespace paso {
31
32 FCT_FluxLimiter::FCT_FluxLimiter(const_TransportProblem_ptr tp)
33 {
34 const dim_t n = tp->transport_matrix->getTotalNumRows();
35 const dim_t blockSize = tp->getBlockSize();
36
37 mpi_info = tp->mpi_info;
38 u_tilde = new double[n];
39 MQ = new double[2*n];
40 R = new double[2*n];
41
42 R_coupler.reset(new Coupler(tp->borrowConnector(), 2*blockSize));
43 u_tilde_coupler.reset(new Coupler(tp->borrowConnector(), blockSize));
44 antidiffusive_fluxes.reset(new SystemMatrix(
45 tp->transport_matrix->type, tp->transport_matrix->pattern,
46 tp->transport_matrix->row_block_size,
47 tp->transport_matrix->col_block_size, true));
48 borrowed_lumped_mass_matrix = tp->lumped_mass_matrix;
49 }
50
51 FCT_FluxLimiter::~FCT_FluxLimiter()
52 {
53 delete[] u_tilde;
54 delete[] MQ;
55 delete[] R;
56 }
57
58 // sets the predictor u_tilde from Mu_tilde by solving M_C * u_tilde = Mu_tilde
59 // and calculates the limiters QP and QN
60 void FCT_FluxLimiter::setU_tilde(const double* Mu_tilde)
61 {
62 const dim_t n = getTotalNumRows();
63 const_SystemMatrixPattern_ptr pattern(getFluxPattern());
64
65 #pragma omp parallel for
66 for (index_t i = 0; i < n; ++i) {
67 const double m = borrowed_lumped_mass_matrix[i];
68 if (m > 0) {
69 u_tilde[i] = Mu_tilde[i]/m;
70 } else {
71 u_tilde[i] = Mu_tilde[i];
72 }
73 }
74
75 // distribute u_tilde
76 u_tilde_coupler->startCollect(u_tilde);
77
78 // calculate
79 // MQ_P[i] = lumped_mass_matrix[i] * max_{j} (\tilde{u}[j]- \tilde{u}[i])
80 // MQ_N[i] = lumped_mass_matrix[i] * min_{j} (\tilde{u}[j]- \tilde{u}[i])
81
82 // first we calculate the min and max of u_tilde in the main block
83 // QP, QN are used to hold the result
84 #pragma omp parallel for
85 for (index_t i = 0; i < n; ++i) {
86 if (borrowed_lumped_mass_matrix[i] > 0) { // no constraint
87 double u_min_i = LARGE_POSITIVE_FLOAT;
88 double u_max_i = -LARGE_POSITIVE_FLOAT;
89 #pragma ivdep
90 for (index_t iptr_ij = pattern->mainPattern->ptr[i];
91 iptr_ij < pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
92 const index_t j=pattern->mainPattern->index[iptr_ij];
93 const double u_j = u_tilde[j];
94 u_min_i = std::min(u_min_i, u_j);
95 u_max_i = std::max(u_max_i, u_j);
96 }
97 MQ[2*i] = u_min_i;
98 MQ[2*i+1] = u_max_i;
99
100 } else {
101 MQ[2*i ] = LARGE_POSITIVE_FLOAT;
102 MQ[2*i+1] = LARGE_POSITIVE_FLOAT;
103 }
104 }
105
106 // complete distribute u_tilde
107 u_tilde_coupler->finishCollect();
108 const double* remote_u_tilde = u_tilde_coupler->borrowRemoteData();
109
110 // now we look at the couple matrix and set the final value for QP, QN
111 #pragma omp parallel for
112 for (index_t i = 0; i < n; ++i) {
113 if (borrowed_lumped_mass_matrix[i] > 0) { // no constraint
114 const double u_i = u_tilde[i];
115 double u_min_i = MQ[2*i];
116 double u_max_i = MQ[2*i+1];
117 #pragma ivdep
118 for (index_t iptr_ij = pattern->col_couplePattern->ptr[i];
119 iptr_ij < pattern->col_couplePattern->ptr[i+1];
120 iptr_ij++) {
121 const index_t j = pattern->col_couplePattern->index[iptr_ij];
122 const double u_j = remote_u_tilde[j];
123 u_min_i = std::min(u_min_i, u_j);
124 u_max_i = std::max(u_max_i, u_j);
125 }
126 MQ[2*i ] = (u_min_i-u_i)*borrowed_lumped_mass_matrix[i];//M_C*Q_min
127 MQ[2*i+1] = (u_max_i-u_i)*borrowed_lumped_mass_matrix[i];//M_C*Q_max
128 }
129 }
130 }
131
132 // starts to update a vector (not given yet) from the antidiffusion fluxes
133 // in flux_limiter->antidiffusive_fluxes (needs u_tilde and Q)
134 void FCT_FluxLimiter::addLimitedFluxes_Start()
135 {
136 const dim_t n = getTotalNumRows();
137 const_SystemMatrixPattern_ptr pattern(getFluxPattern());
138 const double* remote_u_tilde = u_tilde_coupler->borrowRemoteData();
139 SystemMatrix_ptr adf(antidiffusive_fluxes);
140
141 #pragma omp parallel for
142 for (dim_t i = 0; i < n; ++i) {
143 double R_N_i = 1;
144 double R_P_i = 1;
145 if (borrowed_lumped_mass_matrix[i] > 0) { // no constraint
146 const double u_tilde_i = u_tilde[i];
147 double P_P_i = 0.;
148 double P_N_i = 0.;
149 const double MQ_min = MQ[2*i];
150 const double MQ_max = MQ[2*i+1];
151 #pragma ivdep
152 for (index_t iptr_ij = pattern->mainPattern->ptr[i];
153 iptr_ij < pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
154 const index_t j = pattern->mainPattern->index[iptr_ij];
155 if (i != j ) {
156 const double f_ij=adf->mainBlock->val[iptr_ij];
157 const double u_tilde_j=u_tilde[j];
158 /* pre-limiter */
159 if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
160 adf->mainBlock->val[iptr_ij]=0;
161 } else {
162 if (f_ij <=0) {
163 P_N_i+=f_ij;
164 } else {
165 P_P_i+=f_ij;
166 }
167 }
168 }
169 }
170
171 // now the couple matrix
172 #pragma ivdep
173 for (index_t iptr_ij = pattern->col_couplePattern->ptr[i];
174 iptr_ij < pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
175 const index_t j=pattern->col_couplePattern->index[iptr_ij];
176 const double f_ij=adf->col_coupleBlock->val[iptr_ij];
177 const double u_tilde_j=remote_u_tilde[j];
178 // pre-limiter
179 if (f_ij * (u_tilde_j-u_tilde_i) >= 0) {
180 adf->col_coupleBlock->val[iptr_ij]=0;
181 } else {
182 if (f_ij <= 0) {
183 P_N_i+=f_ij;
184 } else {
185 P_P_i+=f_ij;
186 }
187 }
188 }
189 /* finally the R+ and R- are calculated */
190 if (P_N_i<0) R_N_i=std::min(1., MQ_min/P_N_i);
191 if (P_P_i>0) R_P_i=std::min(1., MQ_max/P_P_i);
192 }
193 R[2*i] = R_N_i;
194 R[2*i+1] = R_P_i;
195 }
196
197 // now we kick off the distribution of the R's
198 R_coupler->startCollect(R);
199 }
200
201
202 // completes the exchange of the R factors and adds the weighted
203 // antidiffusion fluxes to the residual b
204 void FCT_FluxLimiter::addLimitedFluxes_Complete(double* b)
205 {
206 const dim_t n = getTotalNumRows();
207 const_SystemMatrixPattern_ptr pattern(getFluxPattern());
208 const_SystemMatrix_ptr adf(antidiffusive_fluxes);
209 const double* remote_R = R_coupler->finishCollect();
210
211 #pragma omp parallel for
212 for (dim_t i = 0; i < n; ++i) {
213 const double R_N_i = R[2*i];
214 const double R_P_i = R[2*i+1];
215 double f_i = b[i];
216
217 #pragma ivdep
218 for (index_t iptr_ij = pattern->mainPattern->ptr[i];
219 iptr_ij < pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
220 const index_t j = pattern->mainPattern->index[iptr_ij];
221 const double f_ij = adf->mainBlock->val[iptr_ij];
222 const double R_P_j = R[2*j+1];
223 const double R_N_j = R[2*j];
224 const double rtmp=(f_ij>=0 ? std::min(R_P_i, R_N_j) : std::min(R_N_i, R_P_j));
225 f_i += f_ij*rtmp;
226 }
227 #pragma ivdep
228 for (index_t iptr_ij=pattern->col_couplePattern->ptr[i];
229 iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
230 const index_t j = pattern->col_couplePattern->index[iptr_ij];
231 const double f_ij = adf->col_coupleBlock->val[iptr_ij];
232 const double R_P_j = remote_R[2*j+1];
233 const double R_N_j = remote_R[2*j];
234 const double rtmp=(f_ij>=0) ? std::min(R_P_i, R_N_j) : std::min(R_N_i, R_P_j);
235 f_i += f_ij*rtmp;
236 }
237 b[i]=f_i;
238 } // end of i loop
239 }
240
241 } // namespace paso
242

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/FluxLimiter.cpp:3531-3826 /branches/lapack2681/paso/src/FluxLimiter.cpp:2682-2741 /branches/pasowrap/paso/src/FluxLimiter.cpp:3661-3674 /branches/py3_attempt2/paso/src/FluxLimiter.cpp:3871-3891 /branches/restext/paso/src/FluxLimiter.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/FluxLimiter.cpp:3669-3791 /branches/stage3.0/paso/src/FluxLimiter.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/FluxLimiter.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/FluxLimiter.cpp:3517-3974 /release/3.0/paso/src/FluxLimiter.cpp:2591-2601 /trunk/paso/src/FluxLimiter.cpp:4257-4344 /trunk/ripley/test/python/paso/src/FluxLimiter.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26