picard

demo.m

-close all; clear all; clc;
+close all; clear all;
 %% PAREMETERS
 x_left = 0.0;
 x_right = 0.1;
-%A = 0.01^2;
-global rho c_p k diffusivity lamma L T_m T_l T_s T_0;
-rho = 4510;
-c_p = 520;
-k   = 16;
+cross_section = 0.01^2;
+run_type = "picard"; %picard, nr
+% 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;
+k   = 16*cross_section;
 diffusivity = k / rho / c_p;
-lamma = 0.388150542167233;
-L = 325;
+l = 325000;
 T_m = 1670;
 T_l = 2000;
 T_s = 1500;
 T_0 = T_s;
+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   = 4;
-
+S   = 20;
+max_nr_iters = 10;
 dt = 1;
-num_els = 1000;
+num_time_steps = 1;
+
+num_els = 400;
 num_nodes = num_els + 1;
-num_time_steps = 2;
 wait_time_plots = 0.5;%seconds
 
+gpoints_ref = [-0.93246 -0.66120 -0.23861 0.23861 0.66120 0.993246];
+gpoints_weights = [0.17132 0.36076 0.6791 0.6791 0.36076 0.17132];
+get_gpoints = @(xl,xr) (xl + xr)/2.0 + (xr-xl)/2.0*gpoints_ref;
 %% PAREMETERS
 Xs = linspace(x_left,x_right,num_nodes);
+x_av = (Xs(2:end)+Xs(1:end-1))/2.0;
 node2dof = [-1, 1:(num_nodes-2), -1];
 % node2dof = 1:(num_nodes);
 mask = find(node2dof>0);
@@ -40,7 +49,7 @@ u_np1 = f0(Xs');
 
 ylims = [T_s-100, T_l+100];
 
-fig = figure;
+fig = figure('Position', [10 10 1810 1210]);
 plot(Xs, u_np1);
 ylim(ylims);
 pause(wait_time_plots);
@@ -54,55 +63,100 @@ for time_iter=1:num_time_steps
   dirichlet_vals = nan([num_nodes,1]);
   dirichlet_vals(1) = T_l;
   dirichlet_vals(end) = T_s;
-  % ASSEMBLY
-  A = sparse(num_dofs,num_dofs);
-  M = sparse(num_dofs,num_dofs);
-  b = zeros([num_dofs,1]);
-  for ielem=1:num_els
-    inodes_global = [ielem, ielem+1];
-    xr = Xs(ielem+1);
-    xl = Xs(ielem);
-    h  = xr - xl;
-    M_loc = h / 6.0 * [2.0, 1.0; 1.0, 2.0];
-    M_loc = M_loc * rho * c_p / dt;
-    A_loc = 1 / h * [1.0, -1.0; -1.0, 1.0];
-    A_loc = A_loc * k;
-    for inode=1:2
-      inodeg = inodes_global(inode);
-      idof = node2dof(inodeg);
-      if (idof<0)
-        continue;
+  for nr_iter=1:max_nr_iters
+    fprintf("About to solve for iter %i of NR, tstep %i \n\n", nr_iter, time_iter);
+    u_i = u_np1;% newton-raphson
+    % ASSEMBLY
+    A  = sparse(num_dofs,num_dofs);
+    M  = sparse(num_dofs,num_dofs);
+    L  = sparse(num_dofs,num_dofs);
+    Lp = sparse(num_dofs,num_dofs);
+    b  = zeros([num_dofs,1]);
+    app_c_p = zeros([num_els,1]);
+    for ielem=1:num_els
+      inodes_global = [ielem, ielem+1];
+      xr = Xs(ielem+1);
+      xl = Xs(ielem);
+      h  = xr - xl;
+      u_av = mean(u_i(inodes_global));
+      app_c_p(ielem ) = rho * l * flp(u_av) + rho * c_p;
+      gpoints = get_gpoints(xl,xr);
+      M_loc  = zeros(2);
+      L_loc  = zeros(2);
+      Lp_loc = zeros(2);
+      A_loc  = zeros(2);
+      b_loc  = zeros([2,1]);
+      phis = get_basis_funcs_p1_d1(xl,xr);
+      grad_phis = get_grad_basis_funcs_p1_d1(xl,xr);
+      for igpoint=1:length(gpoints) 
+        weight = gpoints_weights(igpoint) * h;
+        x = gpoints(igpoint);
+        u_gp = sum([phis{1}(x), phis{2}(x)]'.*u_i(inodes_global));
+        for inode=1:2
+          for jnode=1:2
+            A_loc(inode,jnode) = A_loc(inode,jnode) + weight * k * grad_phis{inode}(x)*grad_phis{jnode}(x);
+            M_loc(inode,jnode) = M_loc(inode,jnode) + weight * rho * c_p * phis{inode}(x)*phis{jnode}(x);
+            L_loc(inode,jnode) = L_loc(inode,jnode) + weight * rho * l * flp(u_gp) * phis{inode}(x)*phis{jnode}(x);
+            Lp_loc(inode,jnode) = Lp_loc(inode,jnode) + weight * rho * l * flpp(u_gp) * u_gp * phis{inode}(x)*phis{jnode}(x);
+          end
+        end
       end
-      for jnode=1:2
-        jnodeg = inodes_global(jnode);
-        jdof = node2dof(jnodeg);
-        if (jdof>0)
-          M(idof,jdof) = M(idof,jdof) + M_loc(inode,jnode);
-          A(idof,jdof) = A(idof,jdof) + A_loc(inode,jnode);
-        else
-          dval = dirichlet_vals(jnodeg);
-          b(idof) = b(idof) - A_loc(inode,jnode)*dval;
-          b(idof) = b(idof) - M_loc(inode,jnode)*u_n(jnodeg);
+      M_loc  = M_loc / dt;
+      L_loc  = L_loc / dt;
+      Lp_loc = Lp_loc / dt;
+      for inode=1:2
+        inodeg = inodes_global(inode);
+        idof = node2dof(inodeg);
+        if (idof<0)
+          continue;
+        end
+        for jnode=1:2
+          jnodeg = inodes_global(jnode);
+          jdof = node2dof(jnodeg);
+          if (jdof>0)
+            M(idof,jdof) = M(idof,jdof) + M_loc(inode,jnode);
+            L(idof,jdof) = L(idof,jdof) + L_loc(inode,jnode);
+            Lp(idof,jdof) = Lp(idof,jdof) + Lp_loc(inode,jnode);
+            A(idof,jdof) = A(idof,jdof) + A_loc(inode,jnode);
+          else
+            dval = dirichlet_vals(jnodeg);
+            b(idof) = b(idof) - M_loc(inode,jnode)*u_n(jnodeg);
+            b(idof) = b(idof) - L_loc(inode,jnode)*u_n(jnodeg);
+            b(idof) = b(idof) - A_loc(inode,jnode)*dval;
+          end
         end
       end
     end
+    % SOLVE
+    b = b + (M+L) * u_n(mask);
+    stiffness_matrix = (M + L) + A;
+    r = b - stiffness_matrix * u_i(mask);
+    r_norm = sqrt(trapz(Xs(mask),r.^2));
+    if run_type=="picard"
+      u_np1(mask) = mldivide(stiffness_matrix, b);%Picard
+    elseif run_type=="nr"
+        stiffness_matrix = stiffness_matrix;
+    else
+      disp("invalid run type");
+      exit();
+    end
+    u_np1(forced_nodes) = dirichlet_vals(forced_nodes);
+    err = sqrt(sum((u_np1-u_exact).^2));
+    fprintf("Residual norm is %g, L2 error is %g\n\n", r_norm, err);
+    clf(fig);
+    plot(Xs, u_np1,'DisplayName','uh');
+    hold on
+    plot(Xs, u_exact,'DisplayName','exact');
+    ylim(ylims);
+    %yyaxis right
+    %plot(x_av, app_c_p,'DisplayName','apparent c_p');
+    legend;
+    pause(wait_time_plots)
   end
-  % SOLVE
-  b = b + M * u_n(mask);
-  stiffness_matrix = M + A;
-  u_np1(mask) = mldivide(stiffness_matrix, b);
-  u_np1(forced_nodes) = dirichlet_vals(forced_nodes);
-  clf(fig);
-  plot(Xs, u_np1,'DisplayName','uh');
-  hold on
-  plot(Xs, u_exact,'DisplayName','exact');
-  ylim(ylims);
-  legend;
-  pause(wait_time_plots)
 end
 
 function [sol] = exact_sol(x,t)
-  global rho c_p k diffusivity lamma L T_m T_l T_s T_0;
+  global rho c_p k diffusivity lamma l T_m T_l T_s T_0;
   interface_pos = locate_interface(t);
   sol = zeros(flip(size(x)));
   sol(x<=interface_pos) = T_l - (T_l - T_m)*erf(x(x<=interface_pos)/(2*sqrt(diffusivity*t)))/erf(lamma);
@@ -110,6 +164,23 @@ function [sol] = exact_sol(x,t)
 end
 
 function [interface_pos] = locate_interface(t)
-  global rho c_p k diffusivity lamma L T_m T_l T_s T_0;
+  global rho c_p k diffusivity lamma l T_m T_l T_s T_0;
   interface_pos = 2*lamma*sqrt(diffusivity*t);
 end
+
+function [phis] = get_basis_funcs_p1_d1(xl,xr)
+  phi_l = @(x) -(x-xr) / (xr-xl);
+  phi_r = @(x) +(x-xl) / (xr-xl);
+  phis = {phi_l, phi_r};
+end
+function [grad_phis] = get_grad_basis_funcs_p1_d1(xl,xr)
+  grad_phi_l = @(x) -1 / (xr-xl);
+  grad_phi_r = @(x) +1 / (xr-xl);
+  grad_phis = {grad_phi_l, grad_phi_r};
+end
+function [y] = transcendental_fun(x)
+  global diffusivity sts stl;
+  y = sqrt(pi)*x - ...
+  (stl / exp(x^2) / erf(x) ) + ...
+  (sts / exp(x^2) / erfc(x) );
+end