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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3022 - (show annotations)
Wed May 5 00:14:45 2010 UTC (9 years, 5 months ago) by gross
File MIME type: text/plain
File size: 4258 byte(s)
typo fixed.

1
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
37
38 err_t Paso_ReactiveSolver_solve(Paso_ReactiveSolver* support, Paso_TransportProblem* fctp, double* u, const double dt, const double* source, Paso_Options* options, Paso_Performance *pp)
39 {
40 const double EXP_LIM_MIN =PASO_RT_EXP_LIM_MIN;
41 const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
42 index_t fail=0;
43 register double d_ii, m_i, x_i, e_i, u_i, F_i;
44 dim_t i;
45 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
46
47 #pragma omp parallel for schedule(static) private(i, d_ii, m_i, x_i, e_i, u_i, F_i)
48 for (i=0;i<n;++i) {
49 d_ii=fctp->reactive_matrix[i];
50 m_i=fctp->lumped_mass_matrix[i];
51 x_i=dt*d_ii/m_i;
52 if (x_i >= EXP_LIM_MAX) {
53 fail=1;
54 } else {
55 F_i=source[i];
56 e_i=exp(x_i);
57 u_i=e_i*u[i];
58 if ( abs(x_i) > EXP_LIM_MIN) {
59 u_i+=F_i/d_ii*(e_i-1.);
60 } else {
61 u_i+=F_i*dt/m_i * (1. + x_i/2); /* second order approximation of ( exp(x_i)-1)/x_i */
62 }
63 u[i]=u_i;
64 }
65 }
66 #ifdef PASO_MPI
67 {
68 index_t fail_loc = fail;
69 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm);
70 }
71 #endif
72 if (fail < 0 ) {
73 return SOLVER_DIVERGENCE;
74 } else {
75 return SOLVER_NO_ERROR;
76 }
77 }
78
79 Paso_ReactiveSolver* Paso_ReactiveSolver_alloc(Paso_TransportProblem* fctp)
80 {
81 Paso_ReactiveSolver* out=NULL;
82 out=MEMALLOC(1,Paso_ReactiveSolver);
83 if (Paso_checkPtr(out)) return NULL;
84 return out;
85 }
86
87 void Paso_ReactiveSolver_free(Paso_ReactiveSolver* in)
88 {
89 if (in!=NULL) {
90 MEMFREE(in);
91 }
92 }
93 double Paso_ReactiveSolver_getSafeTimeStepSize(Paso_TransportProblem* fctp)
94 {
95
96 const double EXP_LIM_MAX =PASO_RT_EXP_LIM_MAX;
97 double dt_max=LARGE_POSITIVE_FLOAT, dt_max_loc;
98 index_t fail_loc, fail;
99 dim_t i;
100 const dim_t n=Paso_SystemMatrix_getTotalNumRows(fctp->transport_matrix);
101 register double d_ii,m_i;
102 if (Paso_noError()) {
103 /*
104 * calculate time step size:
105 */
106 dt_max=LARGE_POSITIVE_FLOAT;
107 fail=0;
108 #pragma omp parallel private(dt_max_loc, fail_loc)
109 {
110 dt_max_loc=LARGE_POSITIVE_FLOAT;
111 fail_loc=0;
112 #pragma omp for schedule(static) private(i,d_ii,m_i)
113 for (i=0;i<n;++i) {
114 d_ii=fctp->reactive_matrix[i];
115 m_i=fctp->lumped_mass_matrix[i];
116 if (m_i > 0) {
117 if (d_ii>0) dt_max_loc=MIN(dt_max_loc, m_i/d_ii);
118 } else {
119 fail_loc=-1;
120 }
121 }
122 #pragma omp critical
123 {
124 dt_max=MIN(dt_max,dt_max_loc);
125 fail=MIN(fail, fail_loc);
126 }
127 }
128 #ifdef PASO_MPI
129 {
130 double rtmp_loc[2], rtmp[2];
131 rtmp_loc[0]=dt_max;
132 rtmp_loc[1]= (double) fail;
133 MPI_Allreduce(rtmp_loc, rtmp, 2, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm);
134 dt_max=rtmp[0];
135 fail = rtmp[1] < 0 ? -1 : 0;
136 }
137 #endif
138 if (fail < 0 ) {
139 Paso_setError(VALUE_ERROR, "Paso_ReactiveSolver_getSafeTimeStepSize: negative mass term detected.");
140 return -1;
141 } else {
142 if (dt_max<LARGE_POSITIVE_FLOAT) dt_max*=0.5*EXP_LIM_MAX;
143 }
144 }
145 return dt_max;
146 }
147 void Paso_ReactiveSolver_initialize(const double dt, Paso_TransportProblem* fctp, Paso_Options* options)
148 {
149 }

  ViewVC Help
Powered by ViewVC 1.1.26