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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3793 - (hide annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 7 months ago) by gross
File MIME type: text/plain
File size: 4420 byte(s)
new implementation of FCT solver with some modifications to the python interface
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 gross 3020
37    
38 gross 3793 err_t Paso_ReactiveSolver_solve(Paso_ReactiveSolver* support, Paso_TransportProblem* fctp, double* u, double* u_old, const double* source, Paso_Options* options, Paso_Performance *pp)
39 gross 2987 {
40 gross 3020 const double EXP_LIM_MIN =PASO_RT_EXP_LIM_MIN;
41     const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
42 gross 3793 const double dt = support->dt;
43 gross 3020 index_t fail=0;
44     register double d_ii, m_i, x_i, e_i, u_i, F_i;
45     dim_t i;
46     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
47    
48 gross 3022 #pragma omp parallel for schedule(static) private(i, d_ii, m_i, x_i, e_i, u_i, F_i)
49 gross 3020 for (i=0;i<n;++i) {
50     d_ii=fctp->reactive_matrix[i];
51     m_i=fctp->lumped_mass_matrix[i];
52     x_i=dt*d_ii/m_i;
53     if (x_i >= EXP_LIM_MAX) {
54     fail=1;
55     } else {
56     F_i=source[i];
57     e_i=exp(x_i);
58 gross 3793 u_i=e_i*u_old[i];
59 gross 3020 if ( abs(x_i) > EXP_LIM_MIN) {
60     u_i+=F_i/d_ii*(e_i-1.);
61     } else {
62     u_i+=F_i*dt/m_i * (1. + x_i/2); /* second order approximation of ( exp(x_i)-1)/x_i */
63     }
64     u[i]=u_i;
65     }
66     }
67 jfenwick 3259 #ifdef ESYS_MPI
68 gross 3020 {
69     index_t fail_loc = fail;
70 gross 3021 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm);
71 gross 3020 }
72     #endif
73     if (fail < 0 ) {
74     return SOLVER_DIVERGENCE;
75     } else {
76     return SOLVER_NO_ERROR;
77     }
78     }
79 gross 2987
80     Paso_ReactiveSolver* Paso_ReactiveSolver_alloc(Paso_TransportProblem* fctp)
81     {
82     Paso_ReactiveSolver* out=NULL;
83     out=MEMALLOC(1,Paso_ReactiveSolver);
84 jfenwick 3259 if (Esys_checkPtr(out)) return NULL;
85 gross 2987 return out;
86     }
87    
88     void Paso_ReactiveSolver_free(Paso_ReactiveSolver* in)
89     {
90     if (in!=NULL) {
91     MEMFREE(in);
92     }
93     }
94     double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
95     {
96 gross 3020
97     const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
98     double dt_max=LARGE_POSITIVE_FLOAT, dt_max_loc;
99     index_t fail_loc, fail;
100     dim_t i;
101     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
102     register double d_ii,m_i;
103 jfenwick 3259 if (Esys_noError()) {
104 gross 3020 /*
105     * calculate time step size:
106     */
107     dt_max=LARGE_POSITIVE_FLOAT;
108     fail=0;
109     #pragma omp parallel private(dt_max_loc, fail_loc)
110     {
111     dt_max_loc=LARGE_POSITIVE_FLOAT;
112     fail_loc=0;
113     #pragma omp for schedule(static) private(i,d_ii,m_i)
114     for (i=0;i<n;++i) {
115     d_ii=fctp->reactive_matrix[i];
116     m_i=fctp->lumped_mass_matrix[i];
117     if (m_i > 0) {
118 gross 3033 if ( d_ii>0 ) dt_max_loc=MIN(dt_max_loc, m_i/d_ii);
119 gross 3020 } else {
120     fail_loc=-1;
121     }
122     }
123     #pragma omp critical
124     {
125 gross 3033 dt_max=MIN(dt_max, dt_max_loc);
126 gross 3020 fail=MIN(fail, fail_loc);
127     }
128     }
129 jfenwick 3259 #ifdef ESYS_MPI
130 gross 3020 {
131     double rtmp_loc[2], rtmp[2];
132     rtmp_loc[0]=dt_max;
133     rtmp_loc[1]= (double) fail;
134     MPI_Allreduce(rtmp_loc, rtmp, 2, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
135     dt_max=rtmp[0];
136     fail = rtmp[1] < 0 ? -1 : 0;
137     }
138     #endif
139     if (fail < 0 ) {
140 jfenwick 3259 Esys_setError(VALUE_ERROR, "Paso_ReactiveSolver_getSafeTimeStepSize: negative mass term detected.");
141 gross 3020 return -1;
142     } else {
143 gross 3033 if (dt_max < LARGE_POSITIVE_FLOAT ) {
144     dt_max*=0.5*EXP_LIM_MAX;
145     } else {
146     dt_max=LARGE_POSITIVE_FLOAT;
147     }
148 gross 3020 }
149     }
150     return dt_max;
151 lgao 3012 }
152 gross 3793 void Paso_ReactiveSolver_initialize(const double dt, Paso_ReactiveSolver* rsolver, Paso_Options* options)
153 gross 3014 {
154 gross 3793 rsolver->dt=dt;
155 gross 3018 }

  ViewVC Help
Powered by ViewVC 1.1.26