nr

demo.m

@@ -3,8 +3,8 @@ close all; clear all;
 x_left = 0.0;
 x_right = 0.1;
 cross_section = 0.01^2;
-run_type = "picard"; %picard, nr
-% cross_section = 1;
+run_type = "Newton-Raphson"; %Picard, Newton-Raphson
+cross_section = 1;
 global rho c_p k diffusivity lamma l T_m T_l T_s T_0 gpoints_ref gpoints_weights stl sts;
 rho = 4510*cross_section;
 c_p = 52;
@@ -19,8 +19,8 @@ stl = c_p*(T_l - T_m) / l;
 sts = c_p*(T_m - T_s) / l;
 lamma = fzero(@transcendental_fun, 0.388150542167233);
 time = 0.0;
-S   = 20;
-max_nr_iters = 10;
+S   = 10;
+max_nr_iters = 20;
 dt = 1;
 num_time_steps = 1;
 
@@ -42,14 +42,16 @@ num_dofs = sum(node2dof>0);
 
 %% INITIAL CONDITION
 f0  = @(x) T_0*ones(size(Xs))';
-fl   = @(tem) 1/2*(tanh(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0))+1.0);
-flp  = @(tem) S/(T_l - T_s)*(1 - (tanh(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0))).^2);
-flpp = @(tem) -4*S^2/(T_l - T_s)^2*(tanh(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0)).*sech(S*2/(T_l - T_s)*(tem - (T_s + T_l)/2.0)).^2);
+Tsl_av = (T_s + T_l)/2.0;
+cte = S*2/(T_l - T_s);
+fl   = @(tem) 1/2*(tanh(cte*(tem - Tsl_av))+1.0);
+flp  = @(tem) cte/2.0*(1 - (tanh(cte*(tem - Tsl_av))).^2);
+flpp = @(tem) -cte^2*(tanh(cte*(tem - Tsl_av)).*sech(cte*(tem - Tsl_av)).^2);
 u_np1 = f0(Xs');
 
 ylims = [T_s-100, T_l+100];
 
-fig = figure('Position', [10 10 1810 1210]);
+fig = figure('Position', [10 10 1210 910]);
 plot(Xs, u_np1);
 ylim(ylims);
 pause(wait_time_plots);
@@ -64,7 +66,7 @@ for time_iter=1:num_time_steps
   dirichlet_vals(1) = T_l;
   dirichlet_vals(end) = T_s;
   for nr_iter=1:max_nr_iters
-    fprintf("About to solve for iter %i of NR, tstep %i \n\n", nr_iter, time_iter);
+    fprintf("About to solve for iter %i of %s, tstep %i \n\n", nr_iter, run_type, time_iter);
     u_i = u_np1;% newton-raphson
     % ASSEMBLY
     A  = sparse(num_dofs,num_dofs);
@@ -132,10 +134,12 @@ for time_iter=1:num_time_steps
     stiffness_matrix = (M + L) + A;
     r = b - stiffness_matrix * u_i(mask);
     r_norm = sqrt(trapz(Xs(mask),r.^2));
-    if run_type=="picard"
+    if run_type=="Picard"
       u_np1(mask) = mldivide(stiffness_matrix, b);%Picard
-    elseif run_type=="nr"
-        stiffness_matrix = stiffness_matrix;
+    elseif run_type=="Newton-Raphson"
+      stiffness_matrix = stiffness_matrix + Lp;
+      delta_T = mldivide(stiffness_matrix, r);
+      u_np1(mask) = u_i(mask) + delta_T;
     else
       disp("invalid run type");
       exit();