Commits

Sergey A. Bozhenkov  committed 451b717

Change initial guesses in radiation contacts

  • Participants
  • Parent commits 5a3c106

Comments (0)

Files changed (2)

File src/BoundaryCondition.cpp

         }
         double getT20(double T10) const
         {
-            double K1 = material1.getHeatConductivity(position, T10);
             double epsilon = getEpsilon();
-            double deriv1 = A1 + B1*T10;
-            double y = (T10*T10)*(T10*T10) + deriv1 * K1/epsilon/SIGMA;
-            if (y<0)
-                y = 0;
-            y = sqrt(sqrt(y));
-
             HelpFunction function(T10, A2, B2, position, material2, epsilon);
             SecantMethod solver(function, eps_rel, eps_abs, maxIter);
-            double T20 = solver.solve(y, y);
-            return T20;
+            return solver.solve(estimateT1(), estimateT2());
         }
-        double estimateT(double T10_old) const
+        double estimateT1() const
         {
             //simply no flux condition
             return -A1/B1;
         }
+        double estimateT2() const
+        {
+            //simply no flux condition
+            return -A2/B2;
+        }
 
         void setEpsRel(double eps) { eps_rel = eps; };
         void setEpsAbs(double eps) { eps_abs = eps; };
     unsigned int numX1 = left.getNumX();
     double dx1 = left.getdx();
     double dx2 = right.getdx();
-    double T_old = left(index_t-1, numX1 - 1);
     // T1m1 = T_left(x-h), T1m2 = T_left(x-2h)
     // T2p1 = T_right(x+h), T2p2 = T_right(x+2h)
     double T1m1 = left(index_t, numX1-2);
     function.setEpsRel(eps_rel/2.);
     function.setEpsAbs(eps_abs/2.);
     function.setMaxIter(maxIter);
-    double T_new = function.estimateT(T_old);
     SecantMethod solver(function, eps_rel, eps_abs, maxIter);
-    double T10 = solver.solve(T_old, T_new);
+    double T10 = solver.solve(function.estimateT1(), function.estimateT2());//T_old, T_new);
     double T20 = function.getT20(T10);
     return make_pair<double, double>(T10, T20);
 }
         }
         double getT20(double T10) const
         {
-            //estimate from the first equation
-            double K1 = material1.getHeatConductivity(position, T10);
             double epsilon = getEpsilon();
-            double eps1 = material1.getEmissivity();
-            double deriv1 = A1 + B1*T10;
-            double y = (T10*T10)*(T10*T10) +
-                deriv1 * K1/epsilon/SIGMA/geomFactor +
-                (1.0 - geomFactor)*eps1/geomFactor/epsilon *
-                    ( (T10*T10)*(T10*T10) - (T0*T0)*(T0*T0) );
-            if (y<0)
-                y = 0;
-            y = sqrt(sqrt(y));
-
             //solve the second equation
             HelpFunction2 function(T0, geomFactor, T10, A2, B2,
                                    position, material2, epsilon);
             SecantMethod solver(function, eps_rel, eps_abs, maxIter);
-            double T20 = solver.solve(y, y);
-            return T20;
+            return solver.solve(estimateT1(), estimateT2());
         }
-        double estimateT(double T10_old) const
+        double estimateT1() const
         {
             //simply no flux condition
             return -A1/B1;
         }
+        double estimateT2() const
+        {
+            //simply no flux condition
+            return -A2/B2;
+        }
 
         void setEpsRel(double eps) { eps_rel = eps; };
         void setEpsAbs(double eps) { eps_abs = eps; };
     unsigned int numX1 = left.getNumX();
     double dx1 = left.getdx();
     double dx2 = right.getdx();
-    double T_old = left(index_t-1, numX1 - 1);
     // T1m1 = T_left(x-h), T1m2 = T_left(x-2h)
     // T2p1 = T_right(x+h), T2p2 = T_right(x+2h)
     double T1m1 = left(index_t, numX1-2);
     function.setEpsAbs(eps_abs/2.);
     function.setMaxIter(maxIter);
     SecantMethod solver(function, eps_rel, eps_abs, maxIter);
-    double T_new = function.estimateT(T_old);
-    double T10 = solver.solve(T_old, T_new);
+    double T10 = solver.solve(function.estimateT1(), function.estimateT2());
     double T20 = function.getT20(T10);
     return make_pair<double, double>(T10, T20);
 

File src/SecantMethod.cpp

 // =========================================================================
 double SecantMethod::solve(double x0, double x1)
 {
+    //lazy call
     if (x1 == x0)
-    {
-        //double f0 = f(x0);
-        //if (fabs(f0) < eps_abs)
-            //return x0;
-        //else
-            x1 = x0*(1.0+100.0*eps_rel) + 20*eps_abs;
-    }
+        x1 = x0*(1.0+100.0*eps_rel) + 20*eps_abs;
 
     unsigned int counter = 0;
     while (counter < maxIter)
         double f0 = f(x0);
         double f1 = f(x1);
         double df = (f1 - f0)/(x1 - x0);
-        if (df == 0)
+        if (df == 0 || x1 == x0)
             throw runtime_error("Zero derivative in secant solver.");
         double x_new = x1 - f1/df;
         double delta = fabs(x_new - x1);