non-linear quasi-static

Author: gtheler

.gitignore

@@ -56,8 +56,9 @@ auto_links.m4
 nbproject
 
 # outputs of "make check"
+*.pvd
 mesh2d.msh
-mesh2d.vtk
+mesh2d*.vtk
 laplace-square.msh
 nafems-t4-result.msh
 nafems-le1*.vtk
@@ -67,51 +68,53 @@ encased_rod*.vtk
 nafems-le*.msh
 spinning-disk-*.msh
 square-f.msh
-square_tmp.vtk
-cube-radial.vtk
+square_tmp*.vtk
+cube-radial*.vtk
 long-bar-*.msh
 map-*.msh
-mesh2d.vtk
+mesh2d*.vtk
 spinning-disk-parallel-*.vtk
 thermal-square*.msh
-thermal-square.vtk
+thermal-square*.vtk
 thermal-two-squares-*.msh
 ray-effect-*-?.msh
 orthotropic-beam-*.msh
-heater-cylinder-inches.vtk
-poisson-square.msh
-reaction-force.vtk
+heater-cylinder-inches*.vtk
+reaction-force*.vtk
 ray-effect-diffusion-*.msh
 bunny*.vtk
 bunny*.msh
 2dpwr-*.vtk
 two-is-2.txt
 circle*.msh
 airfoil-converged.msh
-airfoil-converged.vtk
+airfoil-converged*.vtk
 circle_perimeter*.dat
 wilson-2d*.vtk
-2dpwr-eighth.vtu
-2dpwr-quarter.vtu
-beam-ldef-hex20.vtu
-beam-ldef-hex27.vtu
-beam-ldef-hex8.vtu
-beam-ldef-tet10.vtu
-beam-ldef-tet4.vtu
+2dpwr-eighth*.vtu
+2dpwr-quarter*.vtu
+beam-ldef-hex20*.vtu
+beam-ldef-hex27*.vtu
+beam-ldef-hex8*.vtu
+beam-ldef-tet10*.vtu
+beam-ldef-tet4*.vtu
 beam-orthotropic-hex20.msh
 beam-orthotropic-hex27.msh
 beam-orthotropic-hex8.msh
 beam-orthotropic-tet10.msh
 beam-orthotropic-tet4.msh
 heater-cylinder-inches-results.msh
-heater-cylinder-inches.vtu
-reaction-force.vtu
-two-cubes-isotropic-functions.vtu
-hex8-ldef.vtk
-hex8-single-stretched-linear.vtk
-hex8-single-stretched-neohookean.vtk
-hex8-single-stretched-svk.vtk
+heater-cylinder-inches*.vtu
+reaction-force*.vtu
+two-cubes-isotropic-functions*.vtu
+hex8-ldef*.vtk
+hex8-single-stretched-linear*.vtk
+hex8-single-stretched-neohookean*.vtk
+hex8-single-stretched-svk*.vtk
 spooles.out
 wilson-2d-*.vtk
 bunny.msh
 Stanford_Bunny.stl
+poisson-square.vtu
+hex8-ldef.txt
+beam-orthotropic-*.vtu

autogen.sh

@@ -90,7 +90,8 @@ echo -n "creating src/Makefile.am... "
 
 cd src/pdes
 
-automatic="// automatically generated by autogen.sh on $(date)"
+# automatic="// automatically generated by autogen.sh on $(date)"
+automatic="// automatically generated by autogen.sh"
 
 echo ${automatic} > methods.h
 echo ${automatic} > parse.c

src/feenox.h

@@ -1958,7 +1958,7 @@ struct feenox_t {
     Mat JM;      // jacobian for mass matrix = M'*phi_dot
     Mat Jb;      // jacobian for rhs vector = dq/dT for both volumetric and BCs
     Mat J_snes;  // jacobian for SNES
-    Mat J_tran;  // jacobian for TS
+    Mat J_ts;    // jacobian for TS
 
     PetscScalar *eigenvalue;    // eigenvalue vector
     Vec *eigenvector;           // eivenvectors vector

src/pdes/build.c

@@ -271,7 +271,11 @@ int feenox_problem_build_element_volumetric(element_t *this) {
   feenox_debug_print_gsl_matrix(feenox.fem.Ki, stdout);
   printf("JKi\n");
   feenox_debug_print_gsl_matrix(feenox.fem.JKi, stdout);
-*/  
+  printf("fi\n");
+  feenox_debug_print_gsl_vector(feenox.fem.fi, stdout);
+  printf("bi\n");
+  feenox_debug_print_gsl_vector(feenox.fem.bi, stdout);
+*/
 
   if (feenox.pde.has_stiffness)  {
     petsc_call(MatSetValues(feenox.pde.K, feenox.fem.current_GJ, l, feenox.fem.current_GJ, l, gsl_matrix_ptr(feenox.fem.Ki, 0, 0), ADD_VALUES));

src/pdes/petsc_snes.c

@@ -1,7 +1,7 @@
 /*------------ -------------- -------- --- ----- ---   --       -            -
  *  feenox's non-linear solver using PETSc routines
  *
- *  Copyright (C) 2020--2021 Jeremy Theler
+ *  Copyright (C) 2020--2025 Jeremy Theler
  *
  *  This file is part of FeenoX <https://www.seamplex.com/feenox>.
  *
@@ -169,7 +169,7 @@ PetscErrorCode feenox_snes_residual(SNES snes, Vec phi, Vec r, void *ctx) {
 }
 
 
-PetscErrorCode feenox_snes_jacobian(SNES snes,Vec phi, Mat J_snes, Mat P, void *ctx) {
+PetscErrorCode feenox_snes_jacobian(SNES snes, Vec phi, Mat J_snes, Mat P, void *ctx) {
 
   // J_snes = (K + JK - Jb)_bc
   // TODO: we want SAME_NONZERO_PATTERN!

src/pdes/petsc_ts.c

@@ -37,6 +37,7 @@ int feenox_problem_solve_petsc_transient(void) {
         feenox_call(feenox_problem_solve_petsc_nonlinear());
         petsc_call(SNESDestroy(&feenox.pde.snes));
         feenox.pde.snes = NULL;
+
       }
     } else {
       feenox_function_to_phi(feenox.pde.initial_condition, feenox.pde.phi);
@@ -50,8 +51,10 @@ int feenox_problem_solve_petsc_transient(void) {
       if (feenox.pde.initial_condition != NULL) {
         feenox_call(feenox_problem_build());
       }  
-      petsc_call(MatDuplicate(feenox.pde.has_jacobian_K ? feenox.pde.JK : feenox.pde.K, MAT_COPY_VALUES, &feenox.pde.J_tran));
-      petsc_call(TSSetIJacobian(feenox.pde.ts, feenox.pde.J_tran, feenox.pde.J_tran, feenox_ts_jacobian, NULL));
+
+      feenox_check_alloc(feenox.pde.J_ts = feenox_problem_create_matrix("J_ts"));
+//      petsc_call(MatDuplicate(feenox.pde.has_jacobian_K ? feenox.pde.JK : feenox.pde.K, MAT_COPY_VALUES, &feenox.pde.J_ts));
+      petsc_call(TSSetIJacobian(feenox.pde.ts, feenox.pde.J_ts, feenox.pde.J_ts, feenox_ts_jacobian, NULL));
 
       petsc_call(TSSetProblemType(feenox.pde.ts, (feenox.pde.math_type == math_type_linear) ? TS_LINEAR : TS_NONLINEAR));
 
@@ -118,11 +121,6 @@ int feenox_problem_setup_ts(TS ts) {
 //  petsc_call(TSSetMaxStepRejections(feenox.pde.ts, 10000));
 //  petsc_call(TSSetMaxSNESFailures(feenox.pde.ts, 1000));
 
-  // options overwrite
-  petsc_call(TSSetFromOptions(ts));    
-  // TODO: this guy complains about DM (?)
-//  petsc_call(TSSetUp(ts));    
-  
   SNES snes;
   petsc_call(TSGetSNES(ts, &snes));
   if (snes != NULL) {
@@ -135,65 +133,88 @@ int feenox_problem_setup_ts(TS ts) {
     }
   }
 
+  // options overwrite
+  petsc_call(TSSetFromOptions(ts));    
 
 
   return FEENOX_OK;
 }
 
 PetscErrorCode feenox_ts_residual(TS ts, PetscReal t, Vec phi, Vec phi_dot, Vec r, void *ctx) {
 
-//  static int count = 0;
-  feenox_var_value(feenox_special_var(t)) = t;
-  
+  feenox_special_var_value(t) = t;
+
+  // TODO: store in a global temporary vector
+  Vec phi_bc;
+  petsc_call(VecDuplicate(phi, &phi_bc));
+  petsc_call(VecCopy(phi, phi_bc));
+
+  feenox_call(feenox_problem_dirichlet_eval());  
   if (feenox.pde.math_type == math_type_nonlinear) {
-    feenox_call(feenox_problem_phi_to_solution(phi, 0));
+    feenox_call(feenox_problem_dirichlet_set_phi(phi_bc));
+    feenox_call(feenox_problem_phi_to_solution(phi_bc, 0));
+    if (feenox.pde.phi_bc != NULL) {
+      petsc_call(VecCopy(phi_bc, feenox.pde.phi_bc));
+    }
   }
 
   // TODO: for time-dependent neumann BCs it should not be needed to re-build the whole matrix, just the RHS
   feenox_call(feenox_problem_build());
 
-  // TODO: be smart about this too
-  feenox_call(feenox_problem_dirichlet_eval());  
-
-  // compute the residual R(t,phi,phi_dot) = M*(phi_dot)_dirichlet + K*(phi)_dirichlet - b
+  // compute the residual R(t,phi,phi_dot) = M*(phi_dot)_dirichlet + K*(phi)_dirichlet - b  (linear)
+  //                                       = M*(phi_dot)_dirichlet + f(phi)_dirichlet - b   (non-linear)
 
-  // TODO: store in a global temporary vector
-  Vec tmp;
-  petsc_call(VecDuplicate(phi, &tmp));
-
+  // TODO: start with the K/f term and add the mass later
   // set dirichlet BCs on the time derivative and multiply by M
   if (feenox.pde.M != NULL) {
-    petsc_call(VecCopy(phi_dot, tmp));
-    feenox_call(feenox_problem_dirichlet_set_phi_dot(tmp));
-    petsc_call(MatMult(feenox.pde.M, tmp, r));
+    Vec phi_dot_bc;
+    petsc_call(VecDuplicate(phi, &phi_dot_bc));
+    petsc_call(VecCopy(phi_dot, phi_dot_bc));
+    feenox_call(feenox_problem_dirichlet_set_phi_dot(phi_dot_bc));
+    petsc_call(MatMult(feenox.pde.M, phi_dot_bc, r));
+    petsc_call(VecDestroy(&phi_dot_bc));
   } else  {
     petsc_call(VecZeroEntries(r));
   }
 
-  // set dirichlet BCs on the solution and multiply by K
-  // TODO: only if the problem does not compute an actual residual!
-  petsc_call(VecCopy(phi, tmp));
-  feenox_call(feenox_problem_dirichlet_set_phi(tmp));
-
-  petsc_call(MatMultAdd(feenox.pde.K, tmp, r, r));
-  petsc_call(VecDestroy(&tmp));
+  // set dirichlet BCs on the solution compute the residual
+  if (feenox.pde.has_internal_fluxes) {
+    // if the problem provides an internal flux, use it
+    petsc_call(VecCopy(feenox.pde.f, r));
+  } else {
+    // otherwise make it up from the stiffness and the solution
+    petsc_call(MatMultAdd(feenox.pde.K, phi_bc, r, r));
+  }
+  petsc_call(VecDestroy(&phi_bc));
 
   petsc_call(VecAXPY(r, -1.0, feenox.pde.b));
 
   // set dirichlet bcs on the residual
   feenox_call(feenox_problem_dirichlet_set_r(r, phi));
-  
+
   return FEENOX_OK;
 }
 
-PetscErrorCode feenox_ts_jacobian(TS ts, PetscReal t, Vec phi, Vec phi_dot, PetscReal s, Mat J, Mat P, void *ctx) {
+PetscErrorCode feenox_ts_jacobian(TS ts, PetscReal t, Vec phi, Vec phi_dot, PetscReal s, Mat J_ts, Mat P, void *ctx) {
 
-  // return (K + s*M)_dirichlet
-  petsc_call(MatCopy(feenox.pde.K, J, SAME_NONZERO_PATTERN));
-  if (feenox.pde.M != NULL) {
-    petsc_call(MatAXPY(J, s, feenox.pde.M, SAME_NONZERO_PATTERN));
+  petsc_call(MatAssemblyBegin(feenox.pde.K, MAT_FINAL_ASSEMBLY));
+  petsc_call(MatAssemblyEnd(feenox.pde.K, MAT_FINAL_ASSEMBLY));
+  petsc_call(MatAssemblyBegin(J_ts, MAT_FINAL_ASSEMBLY));
+  petsc_call(MatAssemblyEnd(J_ts, MAT_FINAL_ASSEMBLY));
+  
+  // return (K + JK - Jb + s*M)_bc
+  petsc_call(MatCopy(feenox.pde.K, J_ts, DIFFERENT_NONZERO_PATTERN));
+  if (feenox.pde.has_jacobian_K) {
+    petsc_call(MatAXPY(J_ts, +1.0, feenox.pde.JK, DIFFERENT_NONZERO_PATTERN));
+  }
+  if (feenox.pde.has_jacobian_b) {
+    petsc_call(MatAXPY(J_ts, -1.0, feenox.pde.Jb, DIFFERENT_NONZERO_PATTERN));
+  }
+  
+  if (feenox.pde.has_mass) {
+    petsc_call(MatAXPY(J_ts, s, feenox.pde.M, DIFFERENT_NONZERO_PATTERN));
   }
-  feenox_call(feenox_problem_dirichlet_set_J(J));
+  feenox_call(feenox_problem_dirichlet_set_J(J_ts));
 
   return FEENOX_OK;
 }

src/pdes/thermal/init.c

@@ -135,7 +135,7 @@ int feenox_problem_init_runtime_thermal(void) {
   thermal.q.non_uniform = feenox_depends_on_space(thermal.q.dependency_variables);
   thermal.q.non_linear = feenox_depends_on_function(thermal.q.dependency_functions, feenox.pde.solution[0]);  
 
-  feenox.pde.has_mass = (feenox_var_value(feenox_special_var(end_time)) > 0) && (feenox.pde.transient_type != transient_type_quasistatic);  
+  feenox.pde.has_mass = (feenox_var_value(feenox_special_var(end_time)) > 0) && (feenox.pde.transient_type != transient_type_quasistatic); 
   if (feenox.pde.has_mass) {
 ///pr_thermal+kappa+usage kappa
 ///pr_thermal+kappa+description Thermal diffusivity in units of area per unit of time. 

tests/.gitignore

@@ -3,6 +3,7 @@
 *.trs
 *.last
 *.m
+*.pvd
 tmp.*
 mesh2d.msh
 square-centered.msh
@@ -84,5 +85,6 @@ wilson-2d-*.vt*
 hex8-ldef.dat
 *.frd
 *.sta
-hex8-ldef.vtk
+hex8-ldef.vt*
 spooles.out
+beam*.vt*

tests/beam-ldef-quasi-many-steps.fee

@@ -0,0 +1,18 @@
+PROBLEM mechanical QUASISTATIC
+READ_MESH beam-cantilever-tet4.msh
+
+ldef = 1
+end_time = 1
+dt_0 = 0.01
+
+E = 100000
+nu = 0.3
+
+BC left  fixed
+BC right tz=-25*t
+
+SOLVE_PROBLEM
+
+IF done
+  PRINT w(0.5,-0.25,+0.01)+0.016062
+ENDIF

tests/beam-ldef-quasi-single-step.fee

@@ -0,0 +1,18 @@
+PROBLEM mechanical QUASISTATIC
+READ_MESH beam-cantilever-tet4.msh
+
+ldef = 1
+end_time = 1
+dt_0 = 1
+
+E = 100000
+nu = 0.3
+
+BC left  fixed
+BC right tz=-25*t
+
+SOLVE_PROBLEM
+
+IF done
+  PRINT w(0.5,-0.25,+0.01)+0.016062
+ENDIF