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

Annotation of /branches/symbolic_from_3470/paso/src/ReactiveSolver.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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 caltinay 3815 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 caltinay 3815 const double dt = support->dt;
43 gross 3020 index_t fail=0;
44 caltinay 3868 register double d_ii, x_i, e_i, u_i, F_i;
45 gross 3020 dim_t i;
46     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
47    
48 caltinay 3868 #pragma omp parallel for schedule(static) private(i, d_ii, x_i, e_i, u_i, F_i)
49 gross 3020 for (i=0;i<n;++i) {
50 caltinay 3868 const double m_i=fctp->lumped_mass_matrix[i];
51     if (m_i>0) {
52    
53     d_ii=fctp->reactive_matrix[i];
54     x_i=dt*d_ii/m_i;
55     if (x_i >= EXP_LIM_MAX) {
56     fail=1;
57     } else {
58     F_i=source[i];
59     e_i=exp(x_i);
60     u_i=e_i*u_old[i];
61     if ( abs(x_i) > EXP_LIM_MIN) {
62     u_i+=F_i/d_ii*(e_i-1.);
63     } else {
64     u_i+=F_i*dt/m_i * (1. + x_i/2); /* second order approximation of ( exp(x_i)-1)/x_i */
65     }
66     u[i]=u_i;
67 gross 3020 }
68 caltinay 3868 } else {
69     u[i]=u_old[i] + dt * source[i] ; /* constraints added */
70 gross 3020 }
71     }
72 jfenwick 3259 #ifdef ESYS_MPI
73 gross 3020 {
74     index_t fail_loc = fail;
75 gross 3021 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm);
76 gross 3020 }
77     #endif
78     if (fail < 0 ) {
79     return SOLVER_DIVERGENCE;
80     } else {
81     return SOLVER_NO_ERROR;
82     }
83     }
84 gross 2987
85     Paso_ReactiveSolver* Paso_ReactiveSolver_alloc(Paso_TransportProblem* fctp)
86     {
87     Paso_ReactiveSolver* out=NULL;
88     out=MEMALLOC(1,Paso_ReactiveSolver);
89 jfenwick 3259 if (Esys_checkPtr(out)) return NULL;
90 gross 2987 return out;
91     }
92    
93     void Paso_ReactiveSolver_free(Paso_ReactiveSolver* in)
94     {
95     if (in!=NULL) {
96     MEMFREE(in);
97     }
98     }
99     double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
100     {
101 gross 3020
102     const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
103 caltinay 3868
104 gross 3020 double dt_max=LARGE_POSITIVE_FLOAT, dt_max_loc;
105     dim_t i;
106     const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
107     register double d_ii,m_i;
108 jfenwick 3259 if (Esys_noError()) {
109 gross 3020 /*
110     * calculate time step size:
111     */
112     dt_max=LARGE_POSITIVE_FLOAT;
113 caltinay 3868 #pragma omp parallel private(dt_max_loc)
114 gross 3020 {
115 caltinay 3868 dt_max_loc=LARGE_POSITIVE_FLOAT;
116 gross 3020 #pragma omp for schedule(static) private(i,d_ii,m_i)
117     for (i=0;i<n;++i) {
118     d_ii=fctp->reactive_matrix[i];
119     m_i=fctp->lumped_mass_matrix[i];
120 caltinay 3868 if (m_i > 0) { /* no constraint */
121 gross 3033 if ( d_ii>0 ) dt_max_loc=MIN(dt_max_loc, m_i/d_ii);
122 gross 3020 }
123     }
124     #pragma omp critical
125     {
126 gross 3033 dt_max=MIN(dt_max, dt_max_loc);
127 gross 3020 }
128     }
129 jfenwick 3259 #ifdef ESYS_MPI
130 gross 3020 {
131 caltinay 3868 dt_max_loc=dt_max;
132     MPI_Allreduce(&dt_max_loc, &dt_max, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
133 gross 3020 }
134     #endif
135 caltinay 3868
136     if (dt_max < LARGE_POSITIVE_FLOAT ) {
137     dt_max*=0.5*EXP_LIM_MAX; /* make sure there is no exp overflow */
138     } else {
139 gross 3033 dt_max=LARGE_POSITIVE_FLOAT;
140 caltinay 3868 }
141 gross 3020 }
142     return dt_max;
143 lgao 3012 }
144 caltinay 3815 void Paso_ReactiveSolver_initialize(const double dt, Paso_ReactiveSolver* rsolver, Paso_Options* options)
145 gross 3014 {
146 caltinay 3815 rsolver->dt=dt;
147 gross 3018 }

  ViewVC Help
Powered by ViewVC 1.1.26