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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3012 - (hide annotations)
Wed Apr 28 03:39:22 2010 UTC (10 years, 6 months ago) by lgao
File MIME type: text/plain
File size: 3859 byte(s)
Fix the MPI compilation bug.


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     /* Paso: reactive solver (D is a diagonal matrix)
18     *
19     * - Mv_t=Dv+q v(0)=u
20     *
21     * to return v(dt)
22     *
23     */
24    
25     /**************************************************************/
26    
27     /* Author: l.gross@uq.edu.au */
28    
29     /**************************************************************/
30    
31    
32     #include "ReactiveSolver.h"
33     #include "PasoUtil.h"
34     #include "Solver.h"
35    
36     err_t Paso_ReactiveSolver_solve(Paso_ReactiveSolver* support, Paso_TransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options, Paso_Performance *pp)
37     {
38     return SOLVER_NO_ERROR;
39     }
40    
41     Paso_ReactiveSolver* Paso_ReactiveSolver_alloc(Paso_TransportProblem* fctp)
42     {
43     Paso_ReactiveSolver* out=NULL;
44     out=MEMALLOC(1,Paso_ReactiveSolver);
45     if (Paso_checkPtr(out)) return NULL;
46     return out;
47     }
48    
49     void Paso_ReactiveSolver_free(Paso_ReactiveSolver* in)
50     {
51     if (in!=NULL) {
52     MEMFREE(in);
53     }
54     }
55     double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
56     {
57 gross 3005 double dt_max=LARGE_POSITIVE_FLOAT;
58    
59    
60    
61     double beta_0=0., beta_0_loc, betap, beta=EPSILON;
62     dim_t i, n;
63     int fail=0;
64     double dt_max_loc;
65     register double m_ii,m_i, d_ii;
66     #ifdef PASO_MPI
67     double rtmp[2], rtmp_loc[2];
68     #endif
69     n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
70    
71     #pragma omp parallel private(beta_0_loc)
72     {
73     beta_0_loc=0.;
74     #pragma omp for schedule(static) private(i,m_ii,m_i, d_ii)
75     for (i=0;i<n;++i) {
76     m_i=fctp->lumped_mass_matrix[i];
77     m_ii=fctp->main_diagonal_mass_matrix[i];
78     /* printf(" %d : %e %e -> %e\n", i, m_i, m_ii, m_i/m_ii); */
79     if ( m_ii>0 ) {
80     beta_0_loc=MAX(m_i/m_ii, beta_0_loc);
81     } else {
82     fail=1;
83     }
84     }
85     #pragma omp critical
86     {
87     beta_0=MAX(beta_0,beta_0_loc);
88     }
89     }
90     #ifdef PASO_MPI
91     rtmp_loc[0]=beta_0;
92 lgao 3012 rtmp_loc[1]=1.*((double) fail );
93 gross 3005 MPI_Allreduce(rtmp_loc, rtmp, 2, MPI_DOUBLE, MPI_MAX, fctp->mpi_info->comm);
94     beta_0=rtmp_loc[0];
95     fail= (rtmp_loc[0] > 0) ? 1 : 0;
96     #endif
97     if (fail>0) {
98     Paso_setError(TYPE_ERROR,"Paso_ReactiveSolver_getSafeTimeStepSize: non-positive main diagonal entry in mass matrix.");
99     }
100    
101     if ( Paso_noError()) {
102     if ( ( beta_0 < 2. ) && ( beta_0 >= 1. ) ) {
103     beta=MAX(2*(beta_0-1.)/beta_0, EPSILON);
104     } else {
105     Paso_setError(TYPE_ERROR,"Paso_ReactiveSolver_getSafeTimeStepSize: non-positive main diagonal entry in mass matrix.");
106     }
107     }
108     printf("beta = %e from beta_0 = %e\n",beta, beta_0-1);
109     if ( Paso_noError()) {
110    
111     #pragma omp parallel private(dt_max_loc, betap)
112     {
113     betap=beta+1.;
114     dt_max_loc=LARGE_POSITIVE_FLOAT;
115     #pragma omp for schedule(static) private(i, m_ii, m_i, d_ii)
116     for (i=0;i<n;++i) {
117     d_ii=fctp->reactive_matrix[i];
118     m_i=fctp->lumped_mass_matrix[i];
119     m_ii=fctp->main_diagonal_mass_matrix[i];
120    
121     if (d_ii > 0.) dt_max_loc=MIN(dt_max_loc, (betap* m_ii - m_i)/d_ii );
122     }
123     #pragma omp critical
124     {
125     dt_max=MIN(dt_max,dt_max_loc);
126     }
127     }
128     dt_max=1./(beta* Paso_Transport_getTheta(fctp));
129     #ifdef PASO_MPI
130     dt_max_loc=dt_max;
131     MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
132     #endif
133     }
134     return dt_max;
135 lgao 3012 }

  ViewVC Help
Powered by ViewVC 1.1.26