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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3793 - (show annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 6 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
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, double* u_old, 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 const double dt = support->dt;
43 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 #pragma omp parallel for schedule(static) private(i, d_ii, m_i, x_i, e_i, u_i, F_i)
49 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 u_i=e_i*u_old[i];
59 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 #ifdef ESYS_MPI
68 {
69 index_t fail_loc = fail;
70 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, fctp->mpi_info->comm);
71 }
72 #endif
73 if (fail < 0 ) {
74 return SOLVER_DIVERGENCE;
75 } else {
76 return SOLVER_NO_ERROR;
77 }
78 }
79
80 Paso_ReactiveSolver* Paso_ReactiveSolver_alloc(Paso_TransportProblem* fctp)
81 {
82 Paso_ReactiveSolver* out=NULL;
83 out=MEMALLOC(1,Paso_ReactiveSolver);
84 if (Esys_checkPtr(out)) return NULL;
85 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
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 if (Esys_noError()) {
104 /*
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 if ( d_ii>0 ) dt_max_loc=MIN(dt_max_loc, m_i/d_ii);
119 } else {
120 fail_loc=-1;
121 }
122 }
123 #pragma omp critical
124 {
125 dt_max=MIN(dt_max, dt_max_loc);
126 fail=MIN(fail, fail_loc);
127 }
128 }
129 #ifdef ESYS_MPI
130 {
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 Esys_setError(VALUE_ERROR, "Paso_ReactiveSolver_getSafeTimeStepSize: negative mass term detected.");
141 return -1;
142 } else {
143 if (dt_max < LARGE_POSITIVE_FLOAT ) {
144 dt_max*=0.5*EXP_LIM_MAX;
145 } else {
146 dt_max=LARGE_POSITIVE_FLOAT;
147 }
148 }
149 }
150 return dt_max;
151 }
152 void Paso_ReactiveSolver_initialize(const double dt, Paso_ReactiveSolver* rsolver, Paso_Options* options)
153 {
154 rsolver->dt=dt;
155 }

  ViewVC Help
Powered by ViewVC 1.1.26